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1 . INTRODUCTION 


In order to provide turbulence models useful for computations of the 
flowfields involved in advanced scramjet combustion systems, a number of fea- 
tures of these flowfields must be considered. These combustion systems involve 
supersonic flows with embedded subsonic regions and recirculation zones, and 
appropriate turbulence models for scramjet applications must address each of 
these. The geometry of advanced combustors is often three-dimensional, so that 
the effects of three-dimensionality in the flowfield on the turbulence charac- 
teristics must be taken into account. Moreover, the combustion process in a 
scramjet system is embedded within a highly turbulent flow, so that the effects 
of turbulence on chemical reaction rates must be considered, particularly, in 
the scramjet context, with respect to ignition phenomena. On the other hand, to 
be of maximum utility in scramjet combustor design, the turbulence modeling 
should be as simple and straightforward as is consonant with the requirements of 
overall accuracy. In this application, predictions of mean flowfield structure, 
the effects of heat release, and mean chemical reaction rates are of greatest 
importance: details of the turbulence structure itself can be approximated if 

the approximations introduced do not materially affect the prediction of overall 
mixing rate, chemical reaction rate, and parameters such as the wall skin fric- 
tion distribution and flowfield pressure gradient. Since it can be expected 
that different effects may dominate in different regions of the flow: non- 
isotropy in recirculation regions; compressibility effects in high speed flow 
regions; and turbulence-chemistry interaction effects in regions in which fuel 
ignition is occurring, a modular approach may be the most efficient turbulence 
model overall. In such an approach, each module contains the turbulence model 
elements which best account for the dominant features of each region of the 
flowfield. 

An assessment of turbulence models for scramjet applications was initiated 
in September 1979. During the first year of this work, as outlined in Ref. 1, 
the major effort involved the examination of the multiple dissipation length 
scale (MDLS) turbulence model, since this approach appeared to offer the 


potential for greater generality than existing models in the context of scramjet- 
related flowfields. In addition to this work, other efforts carried out during 
the first year of this program included the definition of a technique for the 
estimation of the initial conditions required by field-equation turbulence 
models (Ref. 1), an examination of the use of a modified dissipation rate 
equation with the basic k-e two-equation turbulence model (Refs. 1,2), and the 
development of a supersonic-flow compressibility correction to the dissipation 
rate equation in the two-equation (or MDLS) approach (Refs. 1, 3). 

Although the results of Ref. 1 indicated that the MDLS model is slightly 
more general than the basic k-e model, the gain is not worth the added cost of 
solving two additional equations. Furthermore, the flowfields considered in 
the analyses reported in Ref. 1, while fundamental to and underlying many of 
the structures found in scramjet flows, did not involve large scale recircula- 
tion regions where the effects of stress nonisotropy become important. Accord- 
ingly, the focus of the second year's work shifted to an assessment of the 
performance of a variety of turbulence models in low-speed and high-speed 
recirculating flows. Thus, in the work described in this report, several 
turbulence models, including the basic two-equation model, the MDLS variant of 
the two-equation approach, and the algebraic stress model (ASM) first reported 
by Rodi (Ref. 4) and further developed by Sindir (Ref. 5) have been applied to 
the prediction of both supersonic jet and supersonic shear layer flows. 

In complex flowfields it is difficult to separate some aspects of the 
turbulence modeling problem from the numerical problems inherent in different 
computational approaches for solving the governing equations describing the 
flow. These aspects include the treatment of wall boundary conditions, the 
algorithms used to generate the finite-difference form of the equations, and 
the algorithms used to provide the finite-difference solution of the governing 
equations for the particular turbulence model chosen. While not an integral 
part of the work carried out under the present program, several such problems 
were encountered and are discussed. 

The basic features of the different turbulence models investigated during 
the current phase of the turbulence modeling assessment program are outlined in 
the next section. These models include the basic two-equation k-c formulation, 
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the MDLS variant of the two-equation approach, the algebraic stress model (ASM), 
and modified versions of the k-e and ASM approaches. Details of the formulation 
of each of these models are included in appendices. Following a description of 
the models, results of the work carried out to assess the performance of the 
models are given in Section 3. The discussion includes observations with res- 
pect to the interaction of the turbulence models with different numerical solu- 
tion techniques. These techniques include a stream-function/vorticity approach 
for the computation of axi symmetric incompressible recirculating flows and a 
time-split MacCormack technique for planar supersonic recirculating flows. Com- 
parisons of the predictions of different turbulence models with data and with 
each other are presented for both planar and axisymmetric low-speed recircula- 
ting flows. Also comparisons of the performance of a two-equation model with 
different eddy viscosity models are presented for a high-speed recirculating 
flow. Following this, conclusions drawn from the work described in this report 
are stated, and work planned to continue the turbulence modeling assessment and 
definition is outlined in Section 4. 
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NOMENCLATURE 



area of control volume surface 

area through which the wall stress is applied 

-uv/k (equation (17)); k p - ((k p - )/ (y p - y E ))y p (equation 

(45)); coefficient in equation (73) 

o 

coefficients in equations (68) to (71) 
exponent in equation (73) 

coefficient in modeled form of <J>. . (equation (11)) 
coefficient in modeled form of <j>. . 9 (equation (14)) 
coefficient in modeled form of <j>. . ' (equation (15)) 

1 J 5 I 

coefficient in modeled form of <J>. . 9 ' (equation (16)) 
coefficient in equation (13) 

2 

skin friction coefficient, t /0.5 p U- 

w in 

coefficient in equation (44) 

~ 2 

coefficient in modeled form of u-u. (equation (25)) 

1 J 

coefficient in equation (34) 
coefficient in equation (72) 
coefficient in equation (34) 
coefficient in equation (20) 
coefficient in equation (19) 
coefficient in equation (34) 
coefficient in equation (34) 
coefficient in equation (24) 

coefficient in modeled production term of e (equation (30)) 
coefficient in modeled destruction term of e (equation (30)) 
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I 


c 


e 


3 


D 


D 


e 



D 


E 


E* 


e 


f 

h 


k 

k p 

k T 

l 

M 


n 

P 

P 


P 


k 



P 

R o 

R 1 

Re 


2 

coefficient in modeled form of vu k (8u.j/9Xj ) (equation (30)) 
material derivative; also jet diameter 

diameter of inlet flowfield, sudden expansion configuration 
term in modeled form of <f>.. 0 (equation (13)) 
diffusive transport of (equation (22)) 
coefficient in logarithmic law of the wall 
exp (<*Re v )/Re v (equation (40)) 

normalized anisotropies of Reynolds stresses, (u^u^. - 2/36 1 - J .k)/k 
(equations (68) to (71)) 

length-scale function (equations (17) and (18)) 

distance of separation between two parallel walls (equation (18)); 
also step height 

turbulent kinetic energy 

production region turbulent kinetic energy 

transfer region turbulent kinetic energy 

characteristic turbulence length scale 

Mach number 

unit normal vector 

static pressure 

production rate of turbulent kinetic energy (equation (23)) 

production rate of turbulent kinetic energy (equation (38)) 

production rate of Reynolds stress tensor (equation (23)) 

instantaneous pressure fluctuations 

radius of pipe before expansion 

radius of pipe after expansion 

Reynolds number 
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Re v 

r 

t 

U 

U c 

U, 


U i 

U in 


jr 

2 


u i 


U i u j 


u i u j u k 


w 

x 

X • 


viscous sublayer Reynolds number (taken as a constant = 20) 
position vector; also polar radius in axi symmetric flows 
time 

streamwise mean velocity 

streamwise mean velocity along flowfield centerline 
time averaged velocity component (i=l,2,3) 
instantaneous velocity component (i=l,2,3) 
maximum inlet velocity 

streamwise mean velocity at dump plane (inlet velocity) 
friction velocity /t /p 

streamwise component of Reynolds normal stress 
instantaneous velocity fluctuations 
time-averaged components of Reynolds stress 
time-averaged triple velocity correlations 
transverse mean velocity 

transverse component of Reynolds normal stress 
lateral component of Reynolds normal stress 
streamwise distance 
displacement vector 
transverse distance 


Greek Symbols 

a weighting factor in equation (8) 

r 2/3pk^/e (equations (68) to (71)) 

6 
£ 
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Kronecker delta 

dissipation rate of turbulent kinetic energy 



£. . 
1J 


K 

K* 

X 


*ij .1 


*1j »2 

♦u.r 

W 

d>. . 

T 1J ,W 


production region dissipation rate 

transfer region dissipation rate 

dissipation rate tensor of Reynolds stresses 

mean dissipation rate for near-wall kinetic energy budget 
(equation (45)) 

the von Karman constant; also wave number 

KC 1/4 

K S 

coefficient in algebraic stress formulation (equations (68) to 
(71)) 

coefficient in equation (45) 
dynamic viscosity of fluid 
turbulent viscosity 
total viscosity, y t + y 
kinematic viscosity of fluid 
density of fluid 
spreading parameter for jets 

Prandtl number for k transport equation (equation (26)) 

Prandtl number for e transport equation (equation (31)) 
shear stress 

fluctuating velocity part of pressure-strain correlation 
(equation (11)) 

mean strain part of pressure-strain correlation (equation (14)) 
near-wall correction to <j> • . q (equation (15)) 
near-wall correction to 4> . • 7 (equation (16)) 

1 J 9 ^ 

near-wall effects in the pressure-strain correlation (equation 

(10)) 
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stream function 


ft vorticity 



W 

w 

X 

y 


value at the node east of P 

value at the cell boundary between E and P 

value at the node north of P 

value at the node where the control volume is centered 

radial direction 

value at the node south of P 

total 

turbulent 

streamwise direction normal stress 
streamwise direction shear stress 

value at the edge of the viscous sublayer; or transverse 
direction 

value at the node west of P 

value at the wall; or lateral direction normal stress 
streamwise direction 
transverse direction 
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2. FLOW EQUATIONS AND TURBULENCE MODELS 


In this section, the theoretical framework used in the turbulent flow 
calculations is outlined. This is accomplished in ten subsections: Section 2.1 

presents the mean flow equations and discusses the consequences of Reynolds time 
averaging; Section 2.2 introduces the Reynolds stress transport equation and 
summarizes the modeling efforts at this level of closure; Sections 2.3 and 2.4 
present alternative closure schemes via the algebraic stress and k-e models, 
.respectively, and compare the relative advantages and drawbacks of each model; 
Sections 2.5 and‘2.6 introduce, respectively, the k and e transport equations 
and the derivation of the modeled forms of these equations for high Reynolds 
number flows; Section 2.7 discusses the modifications introduced to the e 
transport equation; Section 2.8 provides a review of the multiple-dissipation 
length scale approach; Section 2.9 presents a new non-equilibrium wall-function 
treatment for near-wall velocity profiles, turbulent kinetic energy budgets and 
dissipation rates; and Section 2.10 gives the two-dimensional forms of the mean 
flow and turbulence model equations used in the computations. 


2.1 MEAN FLOW EQUATIONS AND REYNOLDS TIME AVERAGING 


The equations of motion in the absence of external force fields take the 
following tensorial form for uniform-density Newtonian fluids: 


at 


^7 D i°j 


1 ap , a u 

p 8x i 3x j L p 


aui au.xi 

a Xj axJJ 


(i) 


where U. is the component of instantaneous velocity in the x^ direction, P is 
the static pressure, and p and y are the fluid density and dynamic viscosity, 
respectively. These three equations coupled with the conservation of mass 
principle 




= 0 


( 2 ) 


form the Navier-Stokes equations that predict the dynamic behavior of turbulent 
as well as laminar flows. However, practical turbulent flows contain a cascade 
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of eddy sizes that represent a wide range of time and length scales. Hence any 
numerical scheme using equations 0) and (2) for turbulent flow simulations 
would require a grid fine enough to resolve even the smallest turbulent motions. 
This, at least with the current generation of digital computers, is not possible 
for the solution of practical problems. Since a description of the instantane- 
ous flowfield is beyond our present computing capability, a good compromise 
would seem to be prediction of a mean flowfield which is either "time" or 
"ensemble" averaged. This approach, first proposed by Osborne Reynolds in the 
late 19th century, is the starting point for most of today's applied turbulence 
work. 

Reynolds suggested a statistical average of the instantaneous velocity 0^ 
with respect to time such that 

1 C T ~ 

U- e lim 2 j J u i dt (3) 

T -> oo -j 


where T is a time interval which is long compared with the largest turbulence 
time scales, but shorter than the period over which the averaged flow quantities 
may vary. Inherent in this definition is the idea that the instantaneous veloc- 
ity 0. (or by the same token any other flow variable) can be divided into a mean, 
U. , and a fluctuating component, u^ , as 

il. = U. + u. (4) 


It follows, on taking the average of each side of (4), that the mean of the 
fluctuating component is identically zero 

i r J 

lim J u. dt = 0 

T « cx -T 1 

Substituting definition (4) and its counterpart for the static pressure 
P(P = P + p) into the instantaneous flow equations (1), and time averaging as 
shown in (3) leads to the mean flow equation known as the Reynolds equation 



3Xj 


1 

P 3x. 


+ 





U -U . 

1 J 


(5) 
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The only difference between the instantaneous flow equations and the Reynolds 

equation is the appearance of the fluctuating velocity correlation tensor, u -u . 

i J 

defined as 


u 


i U j 


1 f J 

lim J u.u. dt 
T -j 1 J 


This term, generally known as the Reynolds or turbulent stress, is actually the 
fluctuating velocity counterpart of the mean velocity local acceleration term 
3U.U./9X.. However, it is traditionally taken to the right-hand side and inter- 

J * J 

preted as a stress rather than an acceleration term. In this form this term 
completely overwhelms its viscous counterpart in most turbulent flows and 
becomes the sole mechanism for diffusive momentum transport. The Reynolds 
equation coupled with the time-averaged form of the continuity equation (2) 


3U i 

3x i 


= 0 


( 6 ) 


now becomes the governing set of equations for turbulent momentum transfer. 

This set of equations, however, is not "closed" due to the appearance of the 
Reynolds stress tensor u-u. which introduces six additional unknowns to raise 

' J ______ 

the total number of variables in the four equations to ten (U., P, u-u-)- Thus 

* 1 J 

additional relationships need to be developed to express u.u. in terms of known 

* J 

or calculable variables. These efforts constitute the subject of turbulence 
model ing. 


2.2 REYNOLDS STRESS TRANSPORT EQUATION 

An equation governing the transport of Reynolds stress, u.u . , can be 

* J 

derived from the Navier-Stokes equations for a fluid of uniform properties and 
no external force fields by multiplying equation 0) by u. and adding to it the 

J 

same equation with suffixes i and j interchanged. Time averaging the resultant 
gives 
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■mu 


3u.u. 3u.u . 

— ]_J. + it — ! — *L - _ , 
3t U k 3x k 


9U, 


3u 4 u 


i u k 3x k + U j U k 3x k 


3U. 


II 


3U, 

2v7 ‘I 

3U . / 

' 3u.j 3u k 
> 3x k 3x j 

+ 3u j 3u k ) 
3x k 3x 1 / 



III 

P / 8u i 

+ 3 U J) 

IV 







_a_ ) rr 

it 1 1 4- 

u .p 

u.p 

^ J 



3x, 


v 


( u • U .Ui + 6 . . 

I i j k jk p 


ik p 


3u .u . 3u. 3u, \ ) 

-& 1 + u j ^7 + u i ^7 ) v 


(7) 


where p is the fluctuating component of the static pressure P and the 6 1 s are 
Kronecker deltas. 


Equation (7) representing the transport of Reynolds stresses along a mean 
streamline can be divided into the five different terms shown above. Some of 
these terms are "exact" in the sense that they are expressed only in terms of 
the stresses and the mean strain rate, and thus do not require modeling (terms 
I and II); others, however, need to be modeled because they either include 
higher-order correlations (term V) or correlations between turbulence quantities 
that are not known or calculable (terms III, IV, and V). The goal of this 
section is to discuss each of these terms separately and provide the modeled 
forms, when needed, to lay the groundwork for the algebraic stress model, which 
is a special case of this transport equation. 


A. Convective Transport, Term I 

This term expresses the rate of change of the Reynolds stresses uTlTT along 
a mean streamline. It is composed of a temporal change, 3u.u ./3t, and a local 

. * J 

acceleration, U k 3 U|Uj/ 3 x k , both of which contain only the stresses and the 
mean flowfield, and thus require no modeling. The temporal term represents a 
time-dependent variation over a period much longer than the interval used in 
the time-averaging and vanishes in a statistically stationary flow. 
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I II 


B. Production, Term II 

This term represents the rate of generation of the Reynolds stress u7u7 
through the interaction of the stress with the mean strain rate. Since turbu- 
lence needs a continuous input of energy to maintain itself, this then can also 
be regarded as the energy drain from the mean flow to turbulence. Usually given 
the symbol p . . . , this term acts like a source in the u-u . transport equation and 

* J * J 

requires no modeling. 

C. Viscous Dissipation, Term III 

This term represents the destruction of the Reynolds stress correlation 
u.u. through viscous action. Being negative definite it behaves like a sink and 

* J 

couterbalances the gains in the level of u^u^ due to the production p^.. 
Traditionally expressed as e. . this term includes correlations between various 

* J 

fluctuating velocity gradients that are neither known or calculable. So far 
only two workable proposals have been made for modeling this term; both of these 
express e.. . in terms of the dissipation rate of turbulent kinetic energy e 
(which is the contracted form (i = j) of e..) given as 

' %J 

v 8u i /3X| < (9u i /8X| ( , + 8u^/3x.). Each proposal makes certain assumptions with 
respect to the character of the dissipating eddies in modeling e..; therefore, 

■ J 

depending on the nature of the flow one or the other of the proposals (or even a 
combination of them) may be desirable. The first proposal, by Daly and Harlow 
(Ref. 6), and Donaldson (Ref. 7), assumes that the dissipating motions have the 
same structure as the energy containing eddies (a hypothesis plausible for low 
Reynolds number flows), and relates e. . to e through the coefficient u -u ./k. 

J ' J 

The second proposal, first suggested by Rotta (Ref. 8), assumes that dissipating 
motions, at least for high Reynolds number flows, are isotropic in character and 
can be expressed as 2/36 • - e. A more general representation is a combination of 

* J 

both proposals such as 

e^j = o,2/3 6 • j e + (1 - a) -jM- e (8) 
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where a is assigned values between 0 and 1 depending on the nature of the flow, 
e.g., a = 1 for high Reynolds number flows. It can also be shown that for high 
Reynolds number flows the term 


/3u i 8u ii+ au^\ 
\ 8x k 3x • 3x k 3x i / 


is negligibly small so that e.. and e reduce to 


3u. 3u 


3u, 


2v 7 r~- and 2v(^—M , respectively. 


3 \ 3x k 


3x, 


D. Pressure-Strain Correlation, Term IV 

This term plays an influential role in the Reynolds stress transport 
equation and has to be modeled with care. The primary function of this correla- 
tion is to change the relative levels of the normal stresses, and to act as a 
source (or a sink) in the shear stress equations. Since it makes no direct 
contribution to the level of turbulence energy but merely redistributes it among 
the normal stresses, it is called a "redistributive" term. 

The first step in modeling this correlation is to establish a Poisson 
equation for the fluctuating pressure. Following Chou (Ref. 9) this is done by 
taking the divergence 3/3x^ of the transport equation for u & : 


2 
m 

and imposing the fluctuating velocity continuity constraint 3u^/3x^ h 0 to 
obtain the Poisson equation 

, „2 d 3u 3U 0 .2 

— = -2 — — — — - — (u u - u u ) 

e 3X 2 3x )l 3x m 3x A 1 m W 

l 

Integrating this equation over a region of fluid, multiplying by 

(3u./3x. + 3u./3x.) and time averaging yields the following expression for the 

* vj J 

pressure-strain correlation; 


3x 


m 


( Vm " 


Vm } 


(9) 


3 | 8u j> .. _ 1 3P 8U £ 

3x £ |3t U m 3x m - - p 3x £ - u m 3x m v 
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(p . . 

1 J .w 

where dV and dS denote, respectively, volume and surface elements about the 
given point, and 3/8n is the normal derivative at the surface. Prime super- 
scripts indicate quantities evaluated at a distance r from the given point. 

This representation of the pressure-strain correlation can be broken down into 

three contributing terms. These are: <J>. . , resulting from purely turbulence 

* J 9 ' 

interactions, <b.. 9 involving interactions between the mean strain rate and 

turbulence, and <f>.. representing the effects of rigid boundaries on both 
1 J »w 

(J>. . and <j>. . Modeling of the pressure-strain correlation thus reduces to the 

« J 5 1 1 J 9 ^ 

modeling of each of these terms, as summarized below. The constants appearing 
in the various models are given in Table 2.1, page 19. 

Modeling of <f>.. : The term <J>. . has long been identified as the only mecha- 

1 J 9 I * J 9 I 

nism in the stress transport equation (7) that could promote a return to a state 
of isotropy. This can be best observed for decaying turbulence where there is 
no appreciable mean strain, and the only term left in equation (7) that could 
equalize the normal stress components and reduce the shear stress is <p-. 

Also, when isotropic two-point correlations are substituted into the definition 
pf -j (equation 10), this term vanishes, providing further evidence of its 
function as a promoter of isotropy. The following linear form for <j>. . , was 

* J 9 * 

first proposed in Ref. 8: 

*1j .1 = " C 1 t K u j " I 6 ij k ) 

where e/k defines a time scale and c^ is a constant to be determined from 
experiment. This simple linear form for <f>.. , is widely accepted and used 

1 J J I 
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despite the fact that, as shown by Bradshaw (Ref- 10), the actual "return to 
isotropy" process is highly non-linear. More sophisticated non-linear forms 
such as Lumley and Khajeh Nouri's proposal (Ref. 11) have been suggested, but 
these have shown no significant improvement over Rotta's proposal. 


Modeling of <j>.. This term in the pressure-strain correlation represents the 
interactions between the mean strain rate and turbulence. Rotta (Ref. 8) 
obtained a simpler form by assuming symmetry properties for the two-point 
correlations and treating the mean velocity gradient 3U^/3x m as constant over 
the region of integration: 


♦id, 2 


1 3U H 

2n dxJ 
m 


2 

3 u m u i 
3r i 3r j 


?!V!i 

3r l 8r i 


d V 

->» 

r 


( 12 ) 


where r 0 and r. are the Cartesian components of the position vector r. A 
workable modeled form for this term was first devised by Launder (Ref. 12) and 
then further refined by Naot, Shavit and Wolfshtein (Ref. 13), and Launder, 
Reece, and Rodi (Ref. 14), who, working independently and using different 
analytical techniques, obtained the same expression for <(>... 2 : + 




ij,2 


(eg + 8) 


n 


(p ij 




(30 c£ - 2) 
55 



(8 c'X - 2) 2 

TT < D ij - I 


(13) 


where 


p. . = 

U 


3U: 


u .u 


+ u.u. 


3u i 


i k 3x k j k 3x k 


3U, 


3U, 


D ij E -\ u i\ WT + u j u k^7 


p = - u.u 


3U. 


i“j 3x . 
J 


+ When equations (11) and (13) are used values of 1.5 and 0.4 are recommended 
for c-| and c^, respectively. 
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A simpler degenerate version that only includes the dominant first term in the 
above equation is also widely used (Gibson and Launder, Ref. 15, Samaraweera, 
Ref. 16). Since all three terms in equation (13) vanish under contraction, the 
redistributive nature of <f>. . 9 is not destroyed by this approximation, given by 

* J 9 L- 


j »2 


(p. . - 4-6 . .p) 

v lj 3 u 1 


(14) 


Modeling of <j>. . : Rigid boundaries affect the flowfield by impeding the 

j J »w 

transfer of energy from the streamwise direction to that normal to the boundary, 
and also reduce the relative magnitude of the shear stress. As shown by Shir 
(Ref. 17) and Gibson and Launder (Ref. 15), these effects (contained in the 
surface integral in equation (10)) can be modeled in the form of near-wall 
correction terms <J>! . , and <j>! . 9 : 


♦ij.i ’ c i k < u k u m V m s ij - 2 Vj Vi ‘ 2 u k u i n k n j> f ( S77> 

near-wall correction to , (15) 

1 J 5 I 

^i j , 2 = c 2 ^km,2 n k n m 6 ij " f ^ik.2 n k n j “ f"- ^ j k , 2 n k n i^ f ^nT?7^ 

near-wall correction to 4> - • 9 (16) 

1 J 9 ^ 


where r^ is the position vector, SL is a characteristic turbulence length scale, 

f is the length scale function, and n is the unit normal vector to the surface. 

These terms diminish with distance from a rigid boundary and become negligible 

at great distances. This behavior is achieved by defining the length-scale 

function f^/n.^) in such a way that it vanishes as 1/ n^r. approaches zero. 

Here f is assumed to be directly proportional to J l/y, with l interpreted as 

3/2 

the dissipation length scale k /e and y is the normal distance from the 
rigid boundary. The constant of proportionality is chosen to render f of 
value unity in near-wall turbulence. Thus, for a single wall, f becomes 


f = 


Jlhs 

</a 3 ' 2 y 


(17) 
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where k is the von Karman constant and a = -uv/k. For flow between two parallel 
walls both of which influence the flow at a given point it is assumed that the 
effects are additive, i.e. 



k 3 / 2 /c 

K 

1 

,3/2 

iTiTy) + (T/h - y) 


(18) 


where h is the distance of separation between the walls. 


E. Diffusive Transport, Term V 


As shown in equation (7) the diffusive transport process includes three 
mechanisms: transport through triple-velocity correlations, transport through 

pressure-fluctuating-velocity correlations, and molecular (viscous) transport. 

At high Reynolds numbers, molecular transport is usually insignificant and is 
not retained in the equation. Diffusive transport through the pressure- 
fluctuating-velocity correlations is also assumed to be negligible as a result of 
Irwin's study of self-preserving jets in adverse pressure gradients (Ref. 18), 
and Hanjalic and Launder 's work on asymmetric plane channel flows. Even though 
Lumley has made some suggestions for modeling this term (Ref. 19), no proven 
models are as yet available. Within these approximations the diffusive trans- 
port term reduces to the triple-velocity correlation which has to be evaluated 
algebraically or from a transport equation. Hanjalic and Launder (Ref. 20) 
arrive at the following form for this correlation after introducing drastic 
simplifications to its transport equation: 


u • u .u, = -c 
i j k s 


u i u £ 


3u j u k 

8x„ 




8u k u l 


3u.u . 


3x, 


Vji 


X £ / 
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Daly and Harlow (Ref. 6) have proposed a considerably simpler representation 
for u.u.u, that only retains the last term in the above equation 

1 J K 

t 


u i u j u k = " c s e u k u £ 


3x„ 


(20) 


Values of 0.11 and 0.20 are recommended for c^ and c s ». 


respectively. 
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Thin shear flow calculations by Launder et al. (Ref. 14) and plane wall flow 
predictions by Reece (Ref. 21) show about equal success with data for both 
versions. 


TABLE 2.1. Recommended Values for Turbulence Model Constants 
Reynolds Stress Transport Equations. 


K 


a 


£ 


C 

C 


8 1 

e 2 


c 


k 


c 


e 3 


c 2 

C 1 

c 2 

C £ 


0.4187 

0.09 

1.00 

< 2 /(c - 

e 2 

1.44 
1.92 
0.22 
0.36 (c 

Go 



)c 


1/2 

U 



) 


1.8 

0.6 


0.5 


0.3 

2.55 


2.3 ALTERNATIVE CLOSURE VIA THE ALGEBRAIC STRESS MODEL 

An orthodox second-order closure would require the solution of a transport 
equation of the form of equation (7) for each of the stress components. Even 
for two-dimensional flows this can be a formidable task, since, in addition to 
the mean flow equations (equations (5) and (6)), five other transport equations 
(for u , v , uv, k, and e) need to be solved. Under certain assumptions, 
however, the stress transport equations can be reduced to a set of 
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algebraic expressions. This idea forms the basis of the algebraic stress model 
to be discussed in this section. 


Following Rodi (Ref. 4) it is noted that the only terms containing gradi- 
ents of Reynolds stresses in equation (7) are those responsible for convective 
and diffusive transport. Therefore, if these gradients can be eliminated, the 
Reynolds stress transport equation reduces to a set of algebraic equations of 
the general form 

Vj = v7 (yy 3 V 3 y k )- 

The convective and diffusive transport of the stresses can be related to the 
turbulent kinetic energy transport rates by noting that 


D u -u • 
1 j 

Dt 


u -u . 
TJ 


Dk 

Dt 


k 2- 
K Dt 


u -u • 
1 j 


If the rate of variation of u,u./k along a streamline is much lower than that of 

1 ^ 

u.u. itself, then 

* J 


D V.i _ U i U i Dk 

Dt k Dt 


( 21 ) 


Similarly, ou.u. = ? d( k) + kc 

1 J K 

where the operator d stands for the net diffusion rate of the quantity in 

parentheses. If the spatial gradient of u.u. is large compared with that of 

^ 

u.u./k, then 

■ J 



u.u, 

D U. Uj = D(k) (22) 

Also by definition (see Section 2.5) 

m: - D(k) = p " e 


Therefore, if equations (21) and (22) are solved for Dk/Dt and d( k) and 
the results are substituted into the above expression, it becomes 
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where the quantity in the brackets is nothing more than the sum of the produc- 
tion, dissipation, and pressure-strain terms in the full Reynolds stress trans- 
port equation, equation (7). Replacing these terms with their modeled forms 
from Section 2.3 (except for the production term, which is exact) produces a set 
of implicit algebraic expressions for the stresses in terms of the mean strain 
rate, turbulent kinetic energy k and its dissipation rate e, and the stresses 
themselves. Hence the final form of the algebraic stress model (ASM) that 
includes near-wall corrections becomes 


u.u. = 7 — - — r (p. • - I- 6. .e + d>. . *i + d>. . „ + d)! . . + d>! . (23) 

i J (p - e) \ l j 3 “ij 6 ' V 1J,1 V 1J,2 Mj,1 MJ,2/ v ' 


where 


9U . 


p = - Ui u k the production rate of kinetic energy 
k 


5 . • = -(u.u 
ij \ i 


3U . 


+ u.u, 


k 3x k j k 3x 
stresses 


d \) ’ 


the production rate of individual Reynolds 


*.. , = -c,e/k (u-u. - 2/3 6.. k), the modeled form of the contribution of 

1 J 5 I I 1 J * J 

fluctuating quantities to the pressure-strain correlation, equation 

OD 

<J>. . 9 = -Cr, (p. - - 2/3 6- -p), the modeled form of the contribution of mean 
strain effects to the pressure-strain correlation, equation (14) 

♦ij.l = c i e/k (“S n k n m S ij - 3/2 Vj Vi - 3/2 v7 n k n j> f (nfrr)* 
wall correction to <J>. . , , equation (15) 

1 J J I 

+ij,2 = c 2 <W n k\ S ij ' 3m ik,2 Vj - 3/2 *jk,2 Vl> f < Wn i r P’ 

wall correction to <b.. 9 , equation (16) 

•J 


Equation (23) is the version of ASM used in the work described in this report. 


The algebraic stress model represents a significant simplification over 
the full Reynolds stress closure and yet it is versatile and is based on a 
plausible derivation. For two-dimensional elliptic flows it requires the 
solution of three implicit algebraic equations for the stresses and two trans- 
port equations for k and e. However, the formulation usually entails a consid- 
erable amount of algebraic manipulation which can be tedious and costly. In 
addition, as discussed in Section 3, special care is needed in incorporating 
these stresses into the mean flow equations to ensure good stability and 
convergence characteristics. 

2.4 ALTERNATIVE CLOSURE VIA THE k-e MODEL 

The k-e model (also known as the two-equation model) developed by Jones 
and Launder (Ref. 22) introduces another degree of simplification to closure of 
the mean flow equations. This model achieves closure by relating the Reynolds 
stresses to the mean strain rate through the Boussinesq approximation. 



The effective (turbulent or eddy) viscosity appearing above, y^, is defined in 
terms of a characteristic length and velocity (an idea apparently borrowed 
from the kinetic theory of gases). If this length is taken as the turbulence 
length scale k 3//2 /e and the velocity as k 1//2 , y^ can be expressed as 

p t = c p p k 2 /e 

where c is a constant of proportionality. Equation (24) is the version of 

y 

the k-e model used in the present calculations. 

Conceptually the use of an effective viscosity for turbulent flows has 
many drawbacks. Firstly, contrary to the requirements of the kinetic theory of 
gases the large energy containing eddies are not rigid bodies which retain their 
identity, and also their "mean free paths" are usually not small compared with 
the flow dimensions. In addition the simple isotropic effective viscosity 
concept breaks down for complex turbulent flows where the shear stress and the 
velocity gradients may have opposite signs, and the effects of the non-equal 
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normal stress components and secondary strain rates may be substantial. 

However, despite all these shortcomings the k-e model has been successful in 
predicting a wide range of thin shear flows and also some complex turbulent 
flows (with minor modifications to the basic form given in equation (24)). The 
simplicity of its formulation and excellent numerical stability characteristics 
have made this model popular in turbulent flow computations. 


2.5 MODELING OF THE k TRANSPORT EQUATION 


Both the k-e viscosity and the ASM require evaluation of the turbulent 
kinetic energy and its dissipation rate to define turbulent time and length 
scales. This and the next section will present the modeled form of the k and e 
transport equations, respectively. 


The turbulent kinetic energy equation can be obtained from the Reynolds 
stress transport equation, equation (7), by setting i = j and dividing by 2 


3U 


+ II _Jt k = - u u 

St u k 3x k i u k 8x k v 3x k 3x k v 3x k 3x i 
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This equation when modeled and reduced to its high Reynolds number form 
becomes 
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3t U k 3x, 
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where within the high Reynolds number approximation e is represented by 

.2.. in i i _i_ _ i — _ 7“ “7 k 3k 


3u . 3u • 


3x 


1 I / - - - 1/ U 1/ 

— 7 ^— ; and u^u k /2 has been taken as c k u^u k — g-— , the form 


(25) 


9u.3u k 

proposed by Daly and Harlow (equation (20)). The terms v v — r — and pu. 6.. 

oX^oX^ i IK 

have been neglected as discussed in Section 2.2, part E. This form of the k 
transport equation, suitable for use with the ASM, can be further simplified for 
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the k-e model by introducing the additional postulate u^u^ = 2/3 6^ k into the 
diffusion term to yield 


9k I, 3k_ 

at u k 3x, 
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t ak 
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(26) 


where a. = 3/2c /c. is the turbulent kinetic energy Prandtl number. 

K jj K 


2.6 MODELING OF THE e TRANSPORT EQUATION 


The dissipation rate of kinetic energy e is evaluated from its transport 
equation formed by multiplying the u i transport equation (equation (9)) by 
2v au./ax, and time averaging: 

1 J 
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2v 3p 9u k 
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(27) 


Before this equation can be solved, however, all the terms appearing on the 
right-hand side have to be modeled in terms of known or calculable variables. 
This is a formidable task in itself since generally no measurements of these 
quantities are available. 

Terms a and b represent, respectively, generation due to vortex stretching 
and secondary generation by the mean flow. Both of these terms are shown by 
Tennekes and Lumley (Ref. 23) to be of negligible importance at high Reynolds 
numbers, and are dropped from the equation. Terms c and d on the other 
hand serve as the primary source (generation due to vortex stretching by 
turbulence) and sink (destruction by viscous action), respectively, for this 
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correlation and become increasingly important at high Reynolds numbers. Most 

4 - 

prefer to model the sum of these terms as suggested by Launder et al . (Ref. 14). 


c 


£ 1 



(28) 


where c and c are two constants evaluated, respectively, by reference to 
near-wall turbulence and decay of grid turbulence. The currently recommended 
values for these constants are given in Table 2.1. 

The diffusion term e is treated by neglecting the pressure-diffusion terms 
and modeling 




as -c — u,u. 


3e 


'£g e k i 3x.j 
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This form was proposed by Hirt (Ref. 24) and also was used by Hanjalic and 
Launder (Ref. 20). The constant c can be expressed in terms of c_ and c_ as 

£, t-i to 

0.36(c - c ). 3 

e 2 1 

When these approximations are introduced into equation (27) the modeled e 
transport equation becomes 
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(30) 


The form of the e transport equation given by equation (30) is suitable for use 


with the ASM. The k-e model version is obtained by approximating u^u^ as 
2/3 6 ik k in the diffusion term to yield 
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+ U 

3t u k 3x 


= c t- p 
k E 1 k 


k + 3x k ( v 9x k + SxJ (31 > 


where is the turbulent kinetic energy dissipation rate Prandtl number equal 
to 


‘t’ 

Each of these terms becomes infinitely large as the Reynolds number approaches 
infinity. However, they have opposite signs and their sum remains finite. 
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This follows directly from the definition of c : 
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2.7 MODIFICATION TO THE e TRANSPORT EQUATION 

Probably the weakest point in the closure by both the ASM and k-e models is 
the determination of the turbulent kinetic energy dissipation rate. Lack of 
measurements make the modeling of the e transport equation a challenging task. 
Since this equation is source dominated, the modeling of the source and sink 
terms is very critical. Equation (28) is an attempt to approximate the sum 
of these quantities by reference to the production rate of turbulent kinetic 
energy p r , the dissipation rate £, and a turbulent time scale k/e. It has 
been argued by Pope (Ref. 25) and Hanjalic and Launder (Ref. 2) that the 
production term in the e equation should be made more sensitive to irrotational 
straining. Both of these schemes propose to incorporate a new term into the e- 
equation that could impart this characteristic. However, the success of these 
methods in some free shear flows did not carry over to more complex recircula- 
ting flow predictions. 

A different approach has been taken by Hanjalic, Launder and Sindir (Ref. 5) 
who suggest replacing p (= - u^Uj (8U^/3xj)) in the e transport equations 
(30) and (31) by 


+ 

The plausibility of this proposal has long been debated since the fine scale 
motions e represents are isotropic in character and should not be directly 
sensitive to the turbulent kinetic energy production rate (which involves large 
energy carrying eddies strained by the mean flow). 


26 


There is no compelling argument to suggest that the effects of mean strain 
should be accountable by a term exactly proportional to the kinetic energy 
generation rate. Equation (32) has been adopted because it displays greater 
sensitivity than the original to streamwise curvature - a characteristic that 
experiments have shown to be desirable. In a straight thin shear flow in 
local equilibrium the proposed form reduces very nearly to that of the "standard" 
e equation. It was hoped that this modification would produce significant 
improvements in recirculating flow predictions where there is substantial stream- 
wise curvature. As discussed in Section 3, this generally proved to be the 
case. The form given in equation (32) was used with both the ASM and the k-e 
models. When equation (32) is used, the turbulence models will be referred to 
as the "modified" ASM and the "modified" k-e, respectively, to differentiate 
from the standard versions of the models. 

2.8 MULTIPLE-SCALE MODELING 

Both the ASM and the k-e approaches are single-point models which adopt 
a single time scale proportional to the turbulence energy turnover time, namely 
k/e. However, it is overly simplistic, at least conceptually, to assume that a 
single time scale can successfully characterize the rates of progress of differ- 
ent turbulent interactions. This realization led Hanjalic, Launder and 
Schiestel (Ref. 3) to develop the multiple-scale approach that is discussed next. 

The key to the new multiple-scale approach is the recognition that while 
the dissipation equation (equation (31)) and the kinetic energy equation 
(equation (26)) both contain production and dissipation terms, these processes 
occur in different spectral regions of the flow. That is, turbulence energy 
production occurs in the larger eddies in the flow, while dissipation phenomena 
involve primarily the smaller scales. Thus, there must be a transfer of energy 
from the larger scales to the smaller, and this transfer can, in certain situa- 
tions, introduce a lag phenomenon, so that turbulence energy production and 
turbulence energy dissipation do not necessarily both increase or decrease in 
the same region of the flow as is implied by equations (26) and (31). 
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To introduce a model in which the evolution of the different scales 
appropriate to the large-eddy production region and the small-eddy dissipation 
region can be accounted for. Launder and co-workers introduced a partitioning 
of the turbulence energy and its dissipation rate, as shown schematically in 
Figure 2.1. In this figure, a partitioning into three regions is shown. For 
wave numbers less than k-j , a production region is defined, characterized by a 
turbulent kinetic energy kp and a dissipation rate £p. This dissipation rate 
controls the transfer of energy through the transfer region < k < k^. For 
wave numbers higher than k^, turbulence energy is dissipated as heat. A 
separate kinetic energy and dissipation rate equation is written for the 
transfer region, characterized by ky and £y, and the production term in the 
kinetic energy equation for the transfer region is equal to the dissipation rate 
Cp in the production region. 

The partitioning of the energy spectrum that is the key feature of the 
multiple-scale model can clearly be carried out as many times as computer 
capacity will allow, but in practice, a partitioning into three regions appears 
to be sufficient (Ref. 3). This requires two sets of transport equations, given 
the assumption (basic to most turbulence modeling) that the mechanisms involved 
i;-, the final dissipation of turbulent kinetic energy into thermal energy are 
capable of accepting all of the energy transferred to them. This assumption is 
the reason that the physical fluid viscosity does not appear in the turbulence 
dissipation rate equations and is supported by the observed Reynolds number 
invariance of fully-turbulent flows. Further, in practice it also is observed 
that the exact point in the wave-number spectrum at which the energy spectrum is 
partitioned does not appear to exert much influence on the results; however, it 
does appear to be influential in initial condition determination (Ref. 1). 

The model equations for the production and transfer region turbulent 
kinetic energy and dissipation rate are similar in form to the k and e transport 
equations discussed in Sections 2.5 and 2.6, respectively. They are: 

Production Region 

a? ♦ u k - «p + 4(( v + ^)k) <33) 
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in which the subscript p refers to the production region and T to the transfer 
region. In this formulation, the turbulent viscosity is given by 


and 


= ( k P + 


k T > 


e p 


(37) 
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This formulation introduces five coefficients, compared to three for the two- 
equation model, but values for several of these coefficients can be inferred 
from two-equation model results and from examination of limiting cases. The 
procedure used to establish the coefficients is described in detail in Ref. 3; 
the results are 



2.2, Cp 

V 1 



c T = 1 .08 , c T = 1.15, c.. = 0.09 

‘l e T l 2 M 


2.9 WALL-FUNCTION TREATMENT 

Most turbulence models including the present versions of the ASM and k-e 
models) are devices for high Reynolds number flows. However, in the vicinity of 
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solid boundaries where the velocities are small, the low Reynolds number effects 
previously neglected become significant and should be accounted for. This can be 
accomplished either by solving the low Reynolds number form of the transport 
equations or by developing wall -functions that introduce these effects into the 
existing high Reynolds number models. Chieng and Launder (Ref. 26) found that 
the first option required vast amounts of computer time due to the slow conver- 
gence characteristics of the low Reynolds number models. On the other hand a 
new wall -function treatment proposed by the same authors was shown to incorpo- 
rate these effects with practically no increase in computing time. An expanded 
version of this treatment is used in the present study. 


A. Near-Wall Velocity Profile and Drag Law 

In flows where the principal source of turbulence energy generation lies 
remote from the walls and the diffusion of energy is towards the surface two key 
approximations can be made in devising wal 1 -functions . Firstly, the dominance 
of the walls on the near-wall length scale can be taken as complete, i.e., 
outside the viscous sublayer (where the flow is viscous but not laminar) the 
turbulence length scale is held to depend, for a limited region near a wall, 
only on the normal distance to the surface. Secondly, the viscous sublayer 

thickness y is assumed to adjust itself according to the external turbulence 

V 1/2 
energy such that the sublayer Reynolds number Re v = y y k y ' /v (where k y is the 

turbulence kinetic energy at the edge of the viscous sublayer) is a universal 

constant equal to 20. Under these conditions Chieng and Launder (Ref. 26) 

propose a new wall -function treatment which is discussed below. 


Figure 2.2a shows a typical near-wall scalar cell (bounded on the west side 
by a wall) on which the following discussions are based. The grid is so 
arranged that node P lies outside the viscous sublayer in the fully turbulent 
region. The shear stress at the wall is estimated by assuming that the mean 
velocity component parallel to the wall varies with height over the fully 
turbulent region as proposed by Launder and Spalding (Ref. 27) 


1/2 

Uk v 

(t w /p) 


&n E* 


yk 


1/2 


(39) 
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with the exception that the kinetic energy is now evaluated at the edge of the 

4 * 

viscous sublayer . Solving this equation for x and evaluating the resultant at 

w 

node P produces the drag law for non-equilibrium turbulent flows 

T w = P % kJ /2 /(£n E* y p kJ /2 /v) (40) 

1/4 1/4 

where k* = <c p and E* = Ec p . E* can also be evaluated in terms of Re y by 
matching the linear velocity profile in the viscous sublayer written as 


Uk 


1/2 


T / O 

w' H 


yk 


1/2 


with the fully turbulent velocity profile (equation (39)) at y v to get 
E* = exp (K*Re v )/Re v . 


B. Near-Wall Turbulent Kinetic Energy Budget and Dissipation Rate 

The near-wall kinetic energy levels are obtained from the solution of the k 
transport equation (equations (25), (26)) modified to reflect these effects. The 
convection and diffusion terms require no changes and are treated in a standard 
way. The production and dissipation terms, however, need to be changed to 
include the near-wall effects. 

The evaluation of the mean production rate of k requires specifying the 
local turbulent shear stress distribution over the near-wall cells. As shown in 
Figure 2.2c a piece-wise continuous shear stress variation is assumed. 


t 


The conventional 


law-of-the-wal 1 


U 


i Ey Vv^ 


valid for flows in local equilibrium (p = e) can be thought„of as a limiting 
case of equation (39) obtained by replacing k with t /pc (the ,1 opal equili- 
brium value of k). This also shows that k* = v kc and l fc* = Ec ' . 
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T t = 


x + (t - t )y/y 


0 i y < y v 

y v i y < y e 


(41) 


This way the stress vanishes in the viscous sublayer, increases abruptly at 
the edge (y = y v ), and then varies linearly over the remainder of the cell. The 
mean production rate of k per unit volume can then be evaluated by integrating 
T t (9U/8y + 9V/9x) over the surface of the cell 


Mean Production 


Rate = y- 
J e y 


/ y< 


/3U . 
T t(ay 



Substituting equation (41) for and equation (39) for U and integrating> 
this expression reduces to its final form 


t (U - U ) 
w v e v' 


\/ T e ~ T w } /, , \ . 3V 

( " y v /y e^ w 3x 


K*pk 



The mean dissipation rate of kinetic energy is evaluated by integrating the 

e distribution over the volume of the cell. In the viscous sublayer the dissi- 

1/2 2 

pation rate is shown to be equal to 2v(3k 7 /3y) by Pope and Whitelaw (Ref. 28), 

and this expression, when coupled with the assumed parabolic variation of 
t 2 

k (= k v (y/y ) ) in the same region becomes 


e = 2vk v /y^ (43) 

In the fully turbulent region, following Spalding (Ref. 29), e is taken to vary 
as 


e 


k 3/2 /c,y 


(44) 


f 

This variation of k corresponds to a linear increase of fluctuating velocities 
with distance from the wall, and also has a zero gradient at the surface. 
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where c 0 is a universal constant given in Table 2.1 as 2.55. The mean dissipa- 
tion rate can now be evaluated by integrating (43) over the viscous Sublayer and 
equation (44) over the fully turbulent region (assuming a linear k variation as 
shown in Figure 2.2b), and averaging the resultant to obtain 


e 



*e c l 




♦ 2a ( k’/ 2 



(45) 


where 


and 



In this equation k p and k £ are the k values at nodes P and E, respectively, and 
k g represents the value at the eastern boundary of the cell as depicted in 
Figure 2.2b. 

The level of k y is obtained by extrapolating the line through k p and k £ to 
y = y v , hence 


k v = k p + ^X (k E- V 


(46) 


The thickness of the viscous sublayer y y and the mean velocity at the location 
U y are then expressed as 

y v - vRe v /kV2 (47) 
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and 


u v = 

It should be noted that y y appears in the above expression for k y 
(equation (46)); hence to obtain a non-iterative solution the following cubic 
equation for (obtained upon substituting for y v from equation (47)) has to be 
solved: 


k 


v 


+ - *y ; 1/2 

- *p 




Alternatively, an iterative scheme that uses the previous iteration level value 
of y v in (46) can be used to evaluate the current k y values. 


2.10 MEAN FLOW AND TURBULENCE MODEL EQUATIONS USED IN THE COMPUTATIONS 

4 * 

The transport and auxiliary equations presented in the previous sections 
reduce, in the case of two-dimensional turbulent flows, to the following 


A. k-e Model 
Conservation of Mass 


3p + 3pU + 1_ 3rpV 
3t 3x r 3r 


x-Momentum 

3pU 3pU 2 1 3rpUV = 

3t 3x r 3r 



(48) 


(49) 


•j* 

These equations are written in a general form to accommodate both Cartesian and 
axi symmetric coordinates. For planar predictions r is set equal to 1 and the 
terms enclosed in boxes are omitted. For axisymmetric flows r equals to the 
polar radius. 
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r-Momentum 


3pV 9pUV 1 9rpV" 
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Turbulence Model 
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and the constants c , c, , c , c and k are defined in Table 2.1. 
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B. Multi-Scale k-e Model 
Conservation of Mass 


3p + ipU + 1 3rpV - n 
3t 3x r 3r 


(54) 


x-Momentum 
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r-Momentum 
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k Transport Equations 
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e Transport Equations 
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where 


2 


+ ip.) + m * / |y. + |v v 

\ 3 x/ \ 3 r/ |\ r / | J \ 3 r 3 x / 


c n = 2.2 


c p - 1.8 - 0.3 [(k p /k T - 1 )/ (k p /k T + 1)] 
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Turbulence Model 


y* = pc - 
K t y 


( k p + k T> i 
2 E P 


( 61 ) 


a. = 1.0 

a = 1.3 

e 

c y = 0.09 
y T = y t + y 


C. Algebraic Stress Model 

The algebraic stress model given below includes wall effects in both the 
the x- and r-directions. Here following Reece (Ref. 21) it is assumed that the 
influence of both sets of walls are simply additive, and there are no cross- 
correlations due to corner effects. Either or both of these wall effects can be 
removed by assigning a value of zero to their corresponding length-scale func- 
tions, f and f . 

A F 


\ a " d f y 


for planar flows. 
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Conservation of Mass 


9£ + ML+ 1 3rpV = n 
3t 3x r 3r 


(62) 


x-Momentum 
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Turbulence Model 


-pu 2 = r[(2 + 4a lx + a lr )(e u + 2/3) + 3a 2x (e u + 2/3) 

- (1 + 2a lx + 2a lr )(e v + 2/3) - 1.5a 2r (e v + 2/3) 

- * 2a 1x - + 2/3 > r 
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where 

a lx = c 2 c 2 ,f x /(1 " c 2^ 

a 2x =C l' tV* 1 - C 2) 
r = 2/3pk 2 /e 

e u = (u 2 - 2/3k)/k 

e w = (w 2 - 2/ 3k ) / k 



a lr = c 2 c 2 ,f r /(1 ' c 2 ) 

a 2r = C 1 ' t V (1 " c 2^ 

A = 0 - c 2 )/(( Cl - 1) + p/e) 
e v - (v 2 - 2/ 3k ) / k 

e uv = ^ /k 

single wall 

parallel walls 
single wall 

parallel walls 
(e.g. coaxial pipes) 


c. , c , c , c , c, , c 9 , Ci 1 and c 9 ' are constants defined in Table 2.1, 
and L and ti are the distances of separation, respectively, between 
the parallel walls in the x- and r-directions. k is the von Karman constant 
(0.4187) and "a" is the near-wall value of -uV/k (generally taken as 0.25). 
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3. ASSESSMENT OF TURBULENCE MODELS FOR SCRAMJET APPLICATIONS 


Three major areas were considered in the assessment of turbulence models 
for scramjet applications carried out under this program. These areas involved 
the development of a compressibility correction technique to extend the accurate 
prediction capability of the thin shear-layer MDLS or k-e models to highly 
supersonic flows; the investigation of the performance capabilities of the k-e, 
MDLS, and ASM models (and the modified k-e and ASM approaches outlined in 
Section 2) in subsonic recirculating flows, both planar and axisymmetric; and 
the development of the k-e and ASM approaches for the prediction of supersonic 
recirculating flows. In the latter two areas problems of turbulence model 
application were encountered that highlight the interaction between turbulence 
modeling and numerical solution techniques, and these problems and their solu- 
tion are described in the subsequent discussions. 

3.1 COMPRESSIBILITY CORRECTION APPROACH 

The basic approach followed in the development of a compressibility correc- 
tion technique involved the introduction of a modification to the term represent- 
ing the dissipation rate production in the dissipation rate expression. Using 
the approximation introduced by Hanjalic, Launder and Sindir and given by 
equation (32), this expression becomes, for thin shear flows 

k P K = p ( c P., c p + C P^ k (frj 

For the MDLS model, Cp = 2.2 and c^ = 0.09. As described in preliminary form in 
Ref. 2, the form of thi correction was taken to be 

c p i = -0.11 + aM b (73) 

for M >_ 1 , where the coefficients a and b are determined from a parametric 
examination of supersonic jet core length and supersonic shear layer growth rate 
predictions. Results of a trial -and-error correlation of supersonic jet core 
length data produced a coefficient fit of the form 
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c pl 1 = -0.11 + 0.0075M, (M> 1) (74) 

that is, a = 0.0075 and b = 1, where M is the local Mach number which varies 
radially and axially throughout the flow. The results of MDLS model predictions 
of jet core lengths are shown in Figure 3.1, as compared to data from a variety 
of sources (Refs. 30-35), and it can be seen from the figure that the increase 
in jet core length as a function of Mach number that occurs for M > 1 is well 
represented. Without the correction given by equation (74), the predicted core 
length trend for M > 1 is equal to that seen in Figure 3.1 for M < 1. Results 
obtained using this correlation with the MDLS model were also compared with the 
supersonic shear layer growth rate data correlation described in Ref. 38, as 
shown on Figure 3.2. While the agreement between the predicted shear layer 
growth rates and those obtained from the experimental data correlation of Ref. 

36 is not as good as for the supersonic jet potential core lengths, the pre- 
dicted results are in considerably better agreement with the overall trend of 
the data using the correction given by equation (74) than are results obtained 
with no correction (Ref. 1), in which there is essentially no change in shear 
layer growth rate with Mach number for M > 1. 

It should be noted that the sensitivity of jet potential core length and 
shear layer growth rate to Mach number exhibited in Figures 3.1 and 3.2 is inde- 
pendent of the observed sensitivity of these parameters to other aspects of the 
flowfield initial conditions. For example, it is well known that potential core 
lengths for subsonic, essentially incompressible round jets are strongly depen- 
dent on the state of flow at the jet exit: laminar, thin boundary layer 

exit conditions produce considerably shorter core lengths than do exit 
conditions which involve thicker, turbulent boundary layers. This phenomenon is 
related to the presence just downstream of the jet exit, in the thin laminar 
boundary layer case, of vortical large scale structures which produce locally 
large mixing rates. These structures do not develop if the jet exit flow 
involves thicker turbulent boundary layers. In the case of the present computa- 
tions (and the experiments to which their results are compared) the initial 
condition in all cases involved a relatively thick turbulent boundary layer at 
the jet nozzle exit. Thus the other aspects of initial condition specification. 
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JET MACH NUMBER, Mj 


FIGURE 3.1. Comparison of Compressibility Correction Results 
With Data for Jet Core Length as a Function of 
Mach Number. 



SPREADING PARAMETER, 
AS DEFINED IN REF. 14 




FIGURE 3.2. Comparison of Compressibility Correction Results 

With Data for Shear Layer Growth Rate as a Function 
of Mach Number. 
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which also affect the determination of velocity potential core lengths, do not 
enter the comparisons just discussed. 

3.2 ASSESSMENT OF MODELS FOR SUBSONIC RECIRCULATING FLOWS 
3.2.1 Planar Flows 

The predictive capabilities of a total of six turbulence models were examined 
for planar backward-facing step flows using the STEP family of elliptic codes 
(Ref. 5). The turbulence models were the k-e model, the multi-dissipation- 
length-scale model, the algebraic stress model (only cross-stream wall effects), 
"modified" k-e model, "modified" algebraic stress model (only cross-stream wall 
effects), and multi -wall algebraic stress model (both cross-stream and stream- 

-L.L 

wise wall effects). These models have been discussed in detail in Chapter 2. 

Before definitive predictions of the test cases were obtained, (i.e., grid- 
independent and fully converged), a preliminary study was conducted to assess 
the relative predictive capability of the models. The mul ti -di ssi pation-length- 
scale (MDLS) model when applied to backward-facing step flow computations failed 
to improve on the standard k-e model predictions. However, the concept of 
multiple time and length scales is physically sound, and is also appealing in 
the sense that it enables separate modeling of turbulent processes proceeding at 
different rates. A more extensive optimization of the coefficients than was 
done by Hanjalic et al . (Ref. 3) may be necessary to extend the predictive capa- 
bility of this model beyond free shear flows. For the work described here the 
MDLS model was abandoned in favor of the k-e model. 

The idea behind the multi -wall ASM formulation was to account for the 
streamwise wall effects on the stress field, and thus try to capture the corner 
eddy that was missing in the previous calculations. Experimental observations 
(Refs. 37-39) indicate a two-eddy structure in all backward-facing step flows as 
shown in Figure 3.3. 


^Influence of walls parallel to the streamwise direction, such as the top and 
bottom walls in backward-facing step geometries. 

^Influence of walls normal to the streamwise direction, such as the rearward 
face of the face. 
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SEPARATED SHEAR LAYER 



FIGURE 3.3. Streamline Pattern for Separating and Reattaching 
Flows in Backward-Facing Step Geometries 



The three dimensional counter-rotating corner eddy which involves lateral as 
well as streamwise and transverse variations extends only a short distance 
(usually less than a step height) downstream of the step, and a large two- 
dimensional primary eddy which involves only streamwise and transverse varia- 
tions occupies the rest of the recirculation zone. Figure 3.3 shows these two 
eddies and the other flow regimes encountered in backward-facing step geometries. 
In this type of flow, probably the one most informative yet simple describing 
parameter is the reattachment length defined as the distance from the step to 

•T* 

the point where the separated shear layer attaches to the surface . This, of 
course, is equal to the combined length of the two eddies. All flow parameters 
in the reverse flow as well as in the adjacent recovery region seem to correlate 
with this quantity, hence an accurate estimate of the reattachment length is 
essential for successful mean flow and turbulence field predictions. Since the 
reattachment length predicted by the ASM generally falls short of the experi- 
mental values by about a step height, it was conjectured that, if this corner 
eddy does grow in magnitude (due to streamwise wall effects), it would "push" 
the reattachment length further downstream and thus bring the computed values 
closer to the measurements. The reattachment length results are given in Table 
3.1 for the Kim, Kline and Johnston 3:1 expansion ratio study (Ref. 39). The 
differences between the multi-wall and standard ASM (i.e., including cross-stream 
wall effects only) appear insignificant. The U-velocity, skin friction and 
peak u7 profiles presented in Figures 3.4 through 3.6 also show no significant 
changes between the two versions of the model. However, contrary to expectations, 
the multi -wall treatment predicts shorter reattachment lengths (by about 2 %) 
and smaller corner eddies. The shortcomings of the multi -wall ASM are probably 
due to the cross correlations between the perpendicular walls that are not 
included in the present version, which adopts Reece's hypothesis (Ref. 21) 
that perpendicular wall effects are simply additive and cross correlations are 
small enough to be neglected. The only other suggestion available in the 
literature, to the best of our knowledge, is Gessner and Eppich's proposal (Ref. 
41) of an "image point" approach. This idea is still in a development stage and 

The point of reattachment is taken as the point where the shear stress 
vanishes. 
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FIGURE 3.4. U-Velocity Profiles at Streamwise Location x/h = 5.333. 

3:1 Area Ratio Planar Expansion, UIN = 18.2 m/s. Data From 
Kim, Kline and Johnston, Ref. 39. 




ABS(uv)/UIN 





is not well tested even for simple shear flows. Therefore, implementation of 
this scheme into recirculating flow computations does not seem appropriate for 
the present. The multi -wall ASM was not pursued further in plane backward- 
facing step flow calculations. 

TABLE 3.1. Reattachment Length Predictions for the Kim, 

Kline and Johnston Study (3:1 Expansion Ratio) 

Reattachment Length 

Models in Step Heights 

ASM 6.16 

Multi -wall ASM 6.02 

Measurements 7 ± 1 

Upon completion of the preliminary study (which eliminated two of the 
original six models) definitive predictions of the test cases were carried out 
using the k-e, "modified" k-e, ASM, and "modified" ASM models for the solution 
domain shown in Figure 3.7. Reattachment length results are given in Table 3.2. 
Predicted and measured profiles of U-velocity and shear stress (uv) are compared 
for three different expansion ratios. These results are now discussed in more 
detail . 

A. U-Velocity Predictions 

Figures 3.8 through 3.10 present, at selected streamwise locations, the 
measured and predicted U-velocity profiles for the 3:1 and 4:1, and the pre- 
dicted profiles for the 9:1 expansion ratios, respectively. A study of these 
results produces two important observations: one is that the relative perfor- 

mance of the models is region-dependent, i.e., best predictions are not 

t ~2 ~~2 ~2 — 

Complete sets of predictions for U, V, P, k, e, u , v , w , uv, p, y. and p/e 

are available for all 3 test cases and 4 models (except for the norma T Reynolds 
stresses with the k-e models) at 40 streamwise stations. To keep the presenta- 
tion manageable profiles for a given variable are presented at representative 
stat-ions and for selected expansion ratios only. However, discussions to follow 
are based on the complete set of results that are kept on file at SAI. Further 
documentation of the results is given in Ref. (5). 
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FIGURE 3.8a. U-Velocity Profiles at Streamwise Location x/h = 2.667. 
3:1 Area Ratio Planar Expansion, UIN = 18.2 m/s. 

Data From Kim, Kline and Johnston, Ref. 39. 
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U-Velocity Profiles at Streamwise Location x/h. = 4.0. 
4:1 Area Ratio Planar Expansion, UIN = 44.2 m/s. 
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FIGURE 3.10b. U-Velocity Profiles at Streamwise Location x/h = 4.0. 

9:1 Area Ratio Planar Expansion, UIN = 44.2 m/s. 












necessarily obtained by the same model in both the reverse flow and recovery 
regions; and secondly, within a given region the established trends do not vary 
with expansion ratio. 

TABLE 3.2 Variation of Reattachment 

4 * 

Length With Expansion Ratio T 

Reattachment Length in Step Heights (h) 

Expansion "Modified" "Modified" 


Ratio 

k-e 

k-e 

ASM 

ASM 

Measurements 

3:1, (40) 

5.36 

5.81 

6.16 

7.77 

7 ± 1 

4:1 ++ , 

4.98 

5.26 

5.40 

6.92 

7.00 ± 0.20 

9:1, (39) 

4.29 

4.56 

4.65 

5.79 

6.25 ± 0.25 


The U-velocity profiles are now discussed starting with the recirculation 
zone. There the "modified" ASM displays the best agreement with data for all 
expansion ratios. The differences between the remaining models are relatively 
small with the ASM, "modified" k-e, and k-e models showing, in that order, 
lesser degrees of agreement. Since they all underestimate the size of the 
reverse flow region as shown in Table 3.2 it is natural that the predicted 
velocities would suffer from "scaling" problems. The maximum mean reverse flow 
velocities, U/U^, are given in Table 3.3 for all models and expansion ratios. 

The magnitude of this velocity seems to decrease slightly with increasing expan- 
sion ratio. Measured values are also reported in that table, however these 
probably do not represent the true maxima since relatively few measuring sta- 
tions were employed in this region (most of them were at least 2 step heights 
apart). Nevertheless, the agreement is surprisingly good. Overall, it was to be 

+ 

Results were obtained with a non-uniform grid of 42 by 42 nodes for a solution 
domain length of 31 step heights ( 1 h before the expansion and 30 downstream of 
the step). Computations were started at x/h = -1 with the measured inlet veloc- 
ity profiles. The inlet kinetic energy and dissipation rates were calculated 
from this profile using Prandtl's mixing length hypothesis. The convergence 
criterion (the maximum acceptable level of the residual sources) was set to 0.001. 

Private communication. Driver, D. M., and Seegmiller, H. L. ; NASA-Ames 
Research Center, STE Branch, 1981. 
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expected that the "modified" ASM would show the best agreement with data, since, 
within the recirculation zone, the U-velocities are closely correlated with the 
reattachment length. This is further supported by the fact that the "modified" 
ASM predicts the lowest spreading rate for the separated shear layer and also 
displays the best agreement with the measured shear stress profiles. 

TABLE 3.3. Peak Mean Reverse Flow Velocities and Locations 


Expansion 

Ratio 


3:1 (a=0°) 


4:1 (a=0°) 


9:1 (a=0°) 


Peak Mean Reverse Location in Step 

Flow Velocity Heights 

U/U x/h 

in 


Model 

Predicted 

+ 

Measured 

Predicted 

Measured 

ASM 

- 0.210 


2.550 


M.ASM 

- 0.230 

- 0.243 

3.210 

5.333 

k-e 

- 0.239 


2.120 


M. k-e 

- 0.245 


2.340 


ASM 

- 0.206 


2.340 


M.ASM 

- 0.229 


3.210 


k-e 

- 0.235 

- 0.213 

1.950 

4.000 

M. k-e 

- 0.240 


2.120 


ASM 

- 0.198 


2.340 


M.ASM 

- 0.229 


2.980 


k-e 

- 0.226 


1.590 


M. k-e 

- 0.235 


1.950 



Downstream of the recirculation zone the predicted recovery of the U- 
velocity profiles is slow compared to the measurements at all expansion ratios. 

^The experimental values given here probably do not represent the true maxima 
since relatively few measuring stations were employed in the recirculation 
region. 
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Here the predictions by the ASM, "modified" k-e and k-e models are hardly dis- 
tinguishable from each other, and in general show fair agreement with data. 

When plotted against the actual downstream distance from the step the "modified" 
ASM seems to suffer from an even slower recovery rate. If the streamwise 
distance is normalized relative to the predicted reattachment points , the 
differences between models diminish. However, the "modified" ASM still displays 
a somewhat different behavior from the other models. Close to the wall it con- 
sistently predicts lower velocities up to an inflection point (at a y/h of about 
0.9) followed by a swift increase (more so than the other models) to the free- 
stream values. This behavior can be partially explained by examining the com- 
puted shear stress profiles. Profiles predicted by the "modified" ASM generally 
reach their peak closer to the wall and drop sooner to the freestream levels 
than the other models. This could explain the more prominent inflective point 
and the faster rise to the freestream levels in the corresponding velocity pre- 
dictions. Close to the wall there is practically no difference in all the pre- 
dicted uv profiles, so the lower velocities of the modified" ASM must be due to 
transport effects (less momentum carry-over due to the larger reverse flow 
region) and to the action of the normal stresses (the magnitude of predicted - 
9(pu 2 )/3x is smaller for the "modified" version). The differences between the 

models in this region diminish with increasing expansion ratio and streamwise 
distance. 

B. Shear Stress Predictions 

Figures 3.11 through 3.13 present at given streamwise locations the 
measured and predicted Ti7 profiles for the 3:1 and 4:1, and the predicted pro- 
files for the 9:1 expansion ratios, respectively. 

An analysis of these figures reveals similar trends in the uv profiles at 
all expansion ratios. Within two step heights downstream of the step, the 
measured uv profiles initially assume small gradually decreasing positive values 
that eventually change sign, and sharply increase in magnitude to a negative 

f 

This takes into account the variation in the starting point of the forward flow 
which may be responsible for the apparent differences between the predicted 
recovery rates. 
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FIGURE 3.11a. uv Profiles at Streamwise Location x/h = 2.333. 

3:1 Area Ratio Planar Expansion, UIN = 18.2 m/s. 

Data From Kim, Kline and Johnston, Ref. 39. 
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FIGURE 3.11b. uv Profiles at Streamwise Location x/h = 5.887. 

3:1 Area Ratio Planar Expansion, UIN = 18.2 m/s. 
Data From Kim, Kline and Johnston, Ref. 39. 
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FIGURE 3.11c. u v Profiles at Streamwise Location x/h = 10.33. 

3:1 Area Ratio Planar Expansion, UIN = 18.2 m/s. 

Data From Kim, Kline and Johnston, Ref. 39. 
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FIGURE 3.12a. 



u7 Profiles at Streamwise Location x/h = 1.25. 

4:1 Area Ratio Planar Expansion, UIN - 44.2 m/s. 

Data From Driver and Seegmiller, Private Communication. 
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FIGURE 3.12b. uv Profiles at Streamwise Location x/h = 4.00. 

4:1 Area Ratio Planar Expansion, UIN = 44.2 m/s. 

Data From Driver and Seegmiller, Private Communication. 
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FIGURE 3.12c. uv Profiles at Streamwise Location x/h = 8.63. 

4:1 Area Ratio Planar Expansion, UIN = 44.2 m/s. 

Data From Driver and Seegmiller, Private Communication. 
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FIGURE 3.13a. uv Profiles at Streamwise Location x/h = 2.0. 

9:1 Area Ratio Planar Expansion, UIN = 44.2 m/s 
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FIGURE 3.13b. uv Profiles at Streamwise Location x/h = 4.00. 

9:1 Area Ratio Planar Expansion, UIN = 44.2 m/s 












peak around a y/h of 0.9. This is followed by a rapid drop to the free stream 
levels. Further downstream the uv values experience a more gradual variation 
with y, increasing monotoni cally (in magnitude) to a negative peak at a y/h of 
about 0.5 before descending to their freestream values. From there, on the same 
qualitative behavior is observed with the profiles flattening out and the 
magnitude of the peaks decreasing with streamwise distance. In all regions 
higher expansion ratios favor a more gradual return to the freestream levels. 

Agreement between predictions and measurements range from satisfactory to 
good with the "modified" ASM generally showing the best agreement followed by 
the ASM, "modified" k-e and k-e models. The differences between the predictions 
diminish downstream of the recirculation zone where all models predict a more 
rapid drop to the freestream levels than the data indicate. This is especially 
true for the "modified" ASM which reaches its peak at lower y/h values and then 
descends to the freestream levels considerably sooner than the other models. 

In all cases the net effect of the modification seems to be a reduction in the 
shear stress levels; maximum peak ov values drop by as much as 52 percent over 
the standard version for the ASM and about 12 percent for the k-e model. This 
decrease also seems to depend on the expansion ratio, with lower expansion 
ratios showing slightly less sensitivity to the modifications in both models. 

The behavior of uv described above brings up several important observations 
for discussion. Firstly, the excellent agreement displayed by the "modified" 

ASM in the recirculation zone is not surprising since this model was also highly 
successful in predicting the size of this region. Similarly, the relative per- 
formance of the other models correlate closely with their reattachment length 
predictions. The effect of the modification on the shear stress computations 
may be explained by studying the expressions used in the models to calculate u7. 
The k-e model employs the Boussinesq approximation (equation (24)) which has a 
transport coefficient Uj defined as c^pk /e. Since one of the effects of the 
modification is to reduce the turbulent kinetic energy, it is reasonable to 
expect that this would also diminish and the shear stress. The ASM, on the 
other hand, uses a more complicated expression (equation (71)) that includes a 
transport coefficient, T (a function of k), and the anisotropic normal stress 
components, e u and e v , to calculate uv. Again the modification is observed to 
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reduce all these quantities leading to lower shear stress levels. The effects 
of the expansion ratio on the uv values can also be explained by reference to 
the turbulent kinetic energy predictions. Since higher expansion ratios tend to 
diminish the levels of k for all models, it is plausible that the u7 values 
which have a strong functional dependence on k would behave similarly. 

C. Conclusions 

All of the models considered successfully predict the qualitative behavior 
of the mean flow and turbulence parameters in all test cases. The relative 
performance of the models however is region-dependent, and hence the best 
predictions are not necessarily obtained by the same model in both the reverse 
flow and recovery zones. The established trends within a given region do not 
vary with expansion ratio or deflection angle. The main effect of these factors 
is to change the size of the recirculation zone and the levels of the turbulent 
kinetic energy and shear stresses. Higher expansion ratios reduce both the size 
of this region and the levels of these quantities. Finally, the modification 
introduced to the models increases the predicted reattachment lengths and 
reduces the turbulent kinetic energy and stress levels in all test cases. In 
the recirculation zone all models tend to underpredict the size of this region 
in varying degrees. Generally, the "modified" ASM computes reattachment lengths 
that compare very favorably with measurements. Predictions by the remaining 
models are usually considerably shorter than the data. Here, the mean velocity 
as well as the turbulence parameters correlate closely with the reattachment 
length. Therefore the best agreement with the measured U- and V-velocities is 
achieved by the "modified" ASM-followed, in order, by the ASM, "modified" k-e, 
and k-e models. The scarcity of reliable turbulence measurements in this region 
makes valid quantitative comparisons a difficult task. Within the scope of 
available data, the agreement between the uv predictions and measurements is 
generally quite good; the "modified" ASM shows the best success -followed, in the 
above order, by the other three models. 

In the recovery region all four models predict too slow a recovery rate. 
Computations by the ASM, "modified" k-e, and k-e models are hardly distinguish- 
able from each other. The "modified" ASM displays an even slower recovery when 
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plotted against the actual distance from the step. When this distance is 
normalized relative to the reattachment point the differences in the predicted 
recovery rates diminish. The agreement between the uv predictions and measure- 
ments is generally satisfactory for all models. 

The following conclusions can be drawn from the assessment of the models 
for plane backward-facing step flows: 

(1) The relative performance of the models in subsonic backward-facing 
step flow computations is region-dependent, i.e., best predictions are not 
necessarily obtained by the same model in both the reverse flow and recovery 
regimes . 

(2) The "modified" ASM produces the best predictions in the reverse flow 
region but computes too slow a recovery rate in the relaxation regime where the 
other models show better agreement with the data. 

(3) A modular approach that uses the "modified" ASM in the reverse flow 
region and the standard version in the recovery regime appears to be the most 
appropriate scheme for predicting subsonic plane backward-facing step flows. 

3.2.2 Axi symmetric Flows 

Before the axi symmetric flow predictions are discussed, it is informative 
to review the past and present efforts in this area. Earlier work which used a 
stream function-vorticity (ip - 9 .) formulation identified three problem areas in 
axi symmetric flow predictions. These were: 

(1) Treatment of the corner point at the expansion plane 

(2) Treatment of the first-derivative term in the tjj-equation 

(3) Treatment of the wall boundary conditions near reattachment 

The corner-point treatment problem involves the proper establishment of 
the boundary conditions at the dump plane, i.e., at the intersection of the 
inlet wall and step. At this point, the dividing streamline that separates the 
recirculation zone downstream of the step from the remainder of the flow has its 
origin. For the proper prediction of recirculation zone length it is essential 
that at this point the dividing streamline have a zero slope, i.e., that it be 
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parallel to the inlet wall contour. Within a stream function-vorticity formula- 
tion, this condition can be directly specified; for primitive-variable formula- 
tions it is specifiable only indirectly, through specification of the gradients 
in the axial and radial velocity components. The second problem involved the 
treatment of the term 1 that appears in the axi symmetric form of the stream 
function equation. As a first derivative term it could not be treated numeri- 
cally in the same fashion as the second-derivative terms which describe the 
remainder of the diffusion of stream function: i.e., in the numerical solution 

of the stream function transport equation 




+ - 
r dr 


-rft 


approximation of the second derivative terms using second-order differences and 
the first derivative term using first-order upwind differencing introduces 
"wiggles" into the solution near the centerline. The problem was overcome by 
noting that 1^- - u (for an incompressible flow) and treating the term in 
question as a "source" term. While this slowed the overall convergence of the 
iterative solution procedure (since the value of the axial mean velocity u is 
obtained for the computation at the last step) this change removed the inaccura- 
cies that resulted from the mixed differencing orders. 

Treatment of the wall boundary condition involves both numerical solution 
technique aspects and turbulence modeling aspects, since within a given overall 
flowfield solution technique it is generally not possible to use a sufficiently 
fine grid to provide the resolution necessary to avoid approximations near 
walls. The standard approach to wall boundary condition formulation is to use 
the law-of-the-wall to define an effective axial velocity at the first grid 
point away from the wall. In general, the logarithmic law of the wall provides 
a valid approximation, but it breaks down near the reattachment point where the 
axial velocity tends toward zero and the logarithmic portion of the boundary 
layer disappears. The solution is to include a treatment, again based on empiri- 
cal formulations, for the laminar sublayer region, and to treat the match point 
between the flowfield solution and the wall boundary as a variable rather than a 
prescribed location. This requires an iteration to be carried out to establish 
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both the friction velocity and the match point, but results in a considerable 
improvement in the prediction of recirculation zone length. 

The effects of the treatment chosen for each of these three problem areas 
in the context of a stream function-vorticity formulation is shown in Figure 
3.14. In all cases the basic two-equation turbulence model, as described in 
Section 2, was used in the flowfield calculation. The flowfield represented is 
an incompressible sudden expansion with two coaxial inlet streams, the inner of 
which is at a velocity level of 1/3 of the outer stream. These data were 
reported by Habib and White! aw (Ref. 41). Substantial differences in the pre- 
dicted centerline velocity results are observed for each of the changes dis- 
cussed in the preceding paragraphs and outlined in Figure 3.14. In each case, 

the same turbulence model has been used, illustrating the effects on the overall 
results of the specification of the numerical aspects of the problem. Since 
different numerical solution procedures involve different boundary condition and 
initial condition treatments, the variation illustrated here can be taken to 
indicate the discrepancies that may be encountered when the turbulence models 
developed using a given solution procedure are applied to analyses which use 
different procedures. 

Much of the computational work described in this section uses a signifi- 
cantly different numerical formulation than the ip-tt code considered in the 

previous paragraphs. This code, which is based on the STEP family, is discussed 

in detail in Ref. 5. It solves two-dimensional (planar or axisymmetric) steady- 
state and time-dependent elliptic partial differential equations through an 
iterative procedure based on an integral control volume analysis with hybrid 

4 * 4 * 4 * 

upwind finite differencing and staggered grids''. The flow and turbulence 

^In evaluating the values of the dependent variables at the cell boundaries one 
of two practices is usually adopted. If it is assumed that the variable in 
question varies linearly between the two nodes bracketing the cell boundary, the 
value of the variable at that boundary is taken as the arithmetical average of 
the nodal values. This procedure is usually referred to as "central differenc- 
ing". The second practice known as "upwind differencing" assumes that the value 
at the cell boundary is equal to that of the upwind node. However, a "hybrid" 
scheme that employs one or the other of these procedures depending on the value 
of the cell Reynolds number may actually be superior to either of these schemes 
used alone. A version of this hybrid upwind finite differencing scheme first 
proposed by Spalding (Ref. 42) is employed in the STEP family. 
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CENTERLINE VELOCITY RATIO, 




FIGURE 3.14. Effects of Code Refinement on Centerline Velocity Profile Prediction. 

Stream Function-Vorticity Code, Two-Equation k-e Turbulence Model. 
Data From Habib and Whitelaw, Ref. 41. 
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equations are solved in primitive variables (U, V, P, k, e, etc.)- Within the 
STEP methodology the treatment of corner cells is straightforward. Figure 3.16 
outlines the treatment used for typical scalar variable, and U- and V-velocity 
corner cells at the expansion plane of a backward-facing step geometry. Again, 
since it is not feasible in the STEP code to provide sufficient near-wall grid 
resolution to avoid approximations, the treatment of wall boundaries near 
reattachment remains a problem in turbulence modeling. In fact, the new non- 
equilibrium wall -functions discussed in Section 2.9 were developed to address 
this particular problem. The basic difference between the two approaches is 
that the STEP family uses the viscous sublayer/fully turbulent region interface 
kinetic energy instead of the friction velocity, while in the stream function- 
vorticity approach the thickness of the viscous sublayer is treated as a variable 
in matching the velocity profiles at the interface between the viscous sublayer 
and the turbulent interior flow. 

Using the STEP code, predictions of the axisymmetric incompressible 2:1 
diameter ratio sudden expansion flowfield described by Chaturvedi (43) were 
obtained with the k-e, "modified" k-e, ASM, and "modified" ASM models. A non- 
uniform 22 by 22 node mesh was used. This mesh is relatively coarse, so that the 
results may not be grid-independent. However, the comparative computations are 
believed to be representative of the capabilities of the different models. To 
provide a broader overall comparison, the results of the predictions carried out 
using the STEP code and the different turbulence models have been compared to 
stream function-vorticity code predictions using the k-e model, and to data 
obtained by Schmotolocha and Phung (Ref. 44) for the same geometry but at a 

r r 

substantially higher inlet Reynolds number (1.4 x 10 compared to 2 x 10 ). 

tt 

For hybrid differencing purposes the STEP family uses a staggered grid in 
which the velocities are evaluated at the boundaries of scalar variable (P, k, e, 
etc.) cells. Hence, separate grids define the locations of the U- and V- 
velocities. A portion of these three grids is shown in Figure 3.15. 

The solution domain is arranged so that the outer surfaces of the boundary 
scalar cells coincide with the physical boundaries of the flow field. The 
advantages of this scheme become apparent when the boundary conditions are 
incorporated into the finite difference formulation. Figure 3.16 shows the 
solution domain used for a typical backward-facing step geometry. 
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FIGURE 3.1 



5. Typical Non-Uniform STEP Grid Portion. 
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Special Treatment for U-Velocity Corner Cell: 



1. Wall shear stress applied only through 
area Ay. 

2. Flux leaves or enters cell from the 
south boundary through area Ag-Ay. 

3. Remaining fluxes calculated in the usual 
way outlined In Ref. 5. 



Special Treatment for V-Velocity Corner Cell: 

1. Wall shear stress applied only through 
area Ay. 

2. Flux leaves or enters cell from the 
west boundary through area A w -Aj. 

3. Remaining fluxes calculated In the usual 
way outlined in Ref. 5. 


CELL *1 


I 1 



Special Treatment for Scalar Variable Corner 

Cell: 

1. For cell *1, flux in the southerly 
direction set to zero. 

2. For cell #2, flux in the westerly 
direction set to zero. 

3. Remaining fluxes calculated in the usual 
way outlined In Ref. 5. 


c) Scalar Variable Corner Cell 
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FIGURE 3.16. Corner Cell Treatment in the STEP Family. 



Reattachment length results are given in Table 3.4. Excluding the "modi- 
fied" ASM which overpredicts the reattachment length by some 40%, the agreement 
between the k-e, "modified" k-e, and ASM computations, and the Chaturvedi measure- 
ments is remarkably good. The k-e model predictions obtained using the ip-Q 
code also show excellent agreement with the data and the STEP code computations. 
The apparent failure of the "modified" ASM is at present difficult to interpret. 
Possible causes are that: the modification to the model may perform poorly in 

f 

axi symmetric flows, or for low expansion ratios , or maybe for both. A series 
of carefully selected test runs to be carried out during continuing work should 
identify the source of the problem. The differences between the remaining three 
models are small. This was anticipated because in low expansion ratio axisym- 
metric geometries, the flow is dominated by the pressure field, and the influ- 
ence of the turbulence models is correspondingly less important. 

TABLE 3.4. Reattachment Length Predictions for the 

Chaturvedi (Ref. 43) Study (2:1 Diameter Ratio) 


Reattachment Length 


Models (STEP Code) in Step Heights 

k-e 7.91 

"modified" k-e 8.07 

ASM 8.08 

"modified" ASM 11.35 

Chaturvedi Data (Ref. 43) - 8.00 

Schmotolocha & Phung Data (Ref. 44) - 9.00 

Stream function-vorticity code predictions 

(k-e model ) 8.00 


Inlet Condition (both computational techniques) 

Reynolds number (based on step height): 2 x 10 5 

Inlet free-stream velocity (U- n ): 31 m/s (see Figure 3.17a) 

Turbulence intensity: 1 x 10“ 5 U? n (assumed uniform) 

+ The lowest expansion ratio tested in planar flows was 3:1. 


85 



Figure 3.17 presents the predicted and measured U-velocity profiles at 
selected streamwise locations. The k-e, "modified" k-e, and ASM models show 
about equal success in predicting the velocity field. The agreement with data 
is generally satisfactory even though all three models tend to underpredict the 
velocities both in the recirculation zone and across the separated shear layer. 
The ASM displays slightly better agreement near the flow centerline, while the 
stream function-vorticity results seem to agree better data in the separated 
shear layer. The "modified" ASM velocity predictions cannot be compared with 
the measurements unless they are properly scaled relative to the predicted reat- 
tachment length. This can be accomplished by defining a new streamwise distance 
as (x^ - x)/h and plotting the predicted velocity profiles in the recirculation 
region relative to this distance, x^ is the computed reattachment length and h 
is the step height. The centerline velocity calculations given in Figure 3.18 
follow the trends established above. The invariance with Reynolds number of 
geometrically similar axi symmetric flows is confirmed by the excellent agreement 
throughout the channel between the radial profiles of velocity measured by 
Chaturvedi (Ref. 43) and by Schmotolocha and Phung (Ref. 44) + - except that 
probe interference effects are evident in the Ref. 44 data near the centerline. 

Profiles of turbulent kinetic energy are compared with Chaturvedi' s 

2 2 

experimental data in Figure 3.19. For these comparisons, the assumptions w = v 

has been made in the experimental data reduction: Chaturvedi presents profiles 

2 2 

of u and v only. At x/R-j = 1.0, both the shape and magnitude of the turbulent 
kinetic energy profiles are predicted extremely well using the stream function- 
vorticity code. The four models tested using the STEP code all successfully sim- 
ulate the shape of the profiles, but they either overpredict the magnitude in the 
case of k-e, "modified" k-e, and ASM models or fall short of the data for the 
"modified" ASM. Further downstream (x/R-j = 2.0 and 3.0), while the magnitude of 
the peak turbulence energy is predicted quite well using the ip-fi formulation and 
reasonably well by the k-e, "modified" k-e and "modified" ASM models using the 
STEP code (the standard ASM still overpredicts the data, progressively greater 
disagreement between the predicted level of turbulence energy and that measured 

Abbott and Kline (Ref. 37) reached the same conclusion for plane backward- 
facing step flows. However, these conclusions only apply to flows that are 
fully turbulent throughout the channel or pipe. 
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FIGURE 3. 1 7d. 


U-Velocity Profiles at x/R-, = 3.00. 2:1 Area Ratio Ax i symmetric 

Expansion, UIN = 31.0 m/s. Data From Chaturvedi , Ref. 43, 
and Schmotolocha & Phung, Ref. 44. 






FIGURE 3.17e. U-Velocity Profiles at x/R-, = 4.00. 2:1 Area Ratio 

Axisymmetric Expansion, UIN = 31.0 m/s. Data From 
Chaturvedi , Ref. 43. 
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by Chaturvedi is evident near the flow centerline. Some disagreement can also 

be seen within the recirculation region. In the latter case at least part of 

2 2 

the disagreement is related to the assumption w = v made in obtaining a turbu- 
lent kinetic energy value from the u 2 and v 2 data presented by Chaturvedi, but 
this is not necessarily the case with respect to the near-centerline turbulent 
kinetic energy levels. In this region the existence of large scale fluctuations 
in an essentially inviscid flow has been observed in other experiments (Ref. 45), 
and these fluctuations, which do not directly contribute to the Reynolds stresses, 
are responsible for the disagreement Between the experiment and the theoretical 
predictions in this region. As mixing proceeds towards the flowfield centerline, 
the contribution of these large-scale fluctuations to the apparent turbulence e 
energy decreases, as is evidenced by the rapid convergence of the experimental 
and predicted turbulent kinetic energy profiles as the flow proceeds downstream. 
Figures 3.19e and 3 . 1 9f . However, some of the energy involved in these large 
scale motions does contribute to the overall flowfield mixing rate, as is shown 
by the velocity profiles in Figures 3.17f and 3.17g: this contribution is not 

adequately accounted for by any of the turbulence models. These results seem to 
indicate that in low expansion ratio axisymmetric flows, the mean velocity and 
turbulence predictions do not change significantly with different turbulence 
models. Simple k-e model predictions appear to be at least as good as the more 
sophisticated ASM computations. 

The effects of ordered large scale fluctuations on the turbulent energy 
profile along the flowfield centerline are shown in a fairly dramatic fashion by 
the centerline turbulent kinetic energy profile comparisons shown in Figure 3.20. 
However, as Figure 3.20 indicates, the overall effect of this phenomenon on the 
centerline mean velocity, and thus the overall mixing rate, is small. Figure 
3.20 can also be taken to show that the overall results indicate an essentially 
inviscid deceleration of the flow up to the end of the recirculation region, 
followed by a turbulent mixing process downstream. The coincidence of the end of 
the region of "inviscid" deceleration and the reattachment point is a feature of 
this particular geometry (a 2:1 radius ratio expansion), since shear layer or jet 
potential core mixing rates are independent of geometry while the reattachment 
length can be correlated directly with the expansion step height. For these pre- 
dictions ASM shows relatively good agreement with the data. The standard and 
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"modified" versions of the k-e model tend to underpredict the magnitude of the 
turbulence kinetic energy but successfully simulate its behavior along the cen- 
terline. Both the "modified" ASM STEP code results and the ip-Q code calculations 
seriously underpredict the data. 

t 

Since this study included only one expansion ratio no definitive 
conclusions (other than the comments made above) can be drawn on the performance 
of the models. A follow-up study covering a range of expansion ratios is needed 
for further evaluation of the models. The discrepancies between the k-e 
model calculations and the STEP code predictions are probably due to the false 
diffusion effects associated with the relatively coarse grid used in the STEP 
code computations. 

3.3 SUPERSONIC RECIRCULATING FLOWS 

Efforts in this area were concentrated on implementing the k-e model 
(Section 2.5) and the non-equilibrium wall-function treatment (Section 2.9) into 
the TWODLE code developed by J. P. Drummond at NASA-Langley (Ref. 46). Reflect- 
ing on the future use of TWODLE it was decided to maintain the built-in turbu- 
lence models in the code and implement the k-e model as an option. This 
required more extensive coding initially but should be beneficial in the sense 
that: 

(1) The original structure of the code (including all subroutines was kept 
intact; therefore, the code can be executed with or without the k-e model. 

(2) Since the model was implemented into a working code (rather than build- 
ing the code around the model) initial debugging was expected to be less time 
consuming. 

(3) This modular approach could also facilitate the incorporation of the 
algebraic stress model which can be introduced as another option. 

These considerations led to the model implementation discussed next. 

4 > 

It should be noted that Chaturvedi data were obtained using hot-wire anemometers 
and relatively primitive techniques. The accuracy of the turbulence data may 
thus be questionable. 

Basically algebraic eddy viscosity specifications. 
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3.3.1 Scheme I 


The initial approach taken was to incorporate the k and e transport equa- 
tions into the MacCormack time-split finite difference technique of TWODLE, and 
solve them simultaneously with the velocity and, if desired, chemistry fields. 
This brought about a basic change in the general form of the TWODLE governing 
equations. A source vector ft was introduced as shown in equation (A1 ) of 
Appendix A to treat the production and dissipation terms associated with the 
k and e equations, or any other source-sink term in the governing equations. 

The new form of the governing equations thus became 



3f^ , 3(a 
3x 3Y 


= H 


where all the terms are defined in Appendix A for the k-c model. 

Within the TWODLE format two ways of incorporating the k and e equations 

into the integration procedure were thought to be possible. One was to solve 

them in the L (At/2) L (At) L (At/2) symmetric operator sequence at both the 
y x y 

intermediate (n + 1) and new (n + 1) time levels with the generation and 
destruction terms calculated only during the L operator step. The solution 

A 

algorithm thus became 

(1) Guess or specify initial fields for all variables. 

(2) Go through the symmetric operator sequence. 

Ly(At/2) / source - sink terms, p(p - e) and pe/k(c p - c c), 

L (At) \ e l e 2 

x' ' { calculated during the operator step L . 

Ly (At/2) 

to obtain the new (n + 1) p, U, V, k and e fields. 

(3) Calculate y. from its definition in terms of c , p, k and e. 

l y 

(4) Go back to step 2 until steady-state is reached. 

The second method was to maintain the practice of solving the k and e transport 
equations in the LJ-^L operator sequence, but calculate the source and sink 
terms at the old (n) time level before the operator sequence was started to 
advance the time. This changed the solution algorithm to 
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(1) Guess or specify initial fields for all variables. 

(2) Calculate y. from its definition in terms of c , p, k and e. 

u y 

(3) Calculate source and sink terms for the k and e equations, 

p(p - e) and pe/k(c p - c c), respectively. 

e l e 2 

(4) Go through the symmetric operator sequence L (At/2) L (At) L (At/2) to 

y * y 

obtain the new (n + 1) p, U, V, k and e fields. 

(5) Go back to step 2 until steady-state is reached. 

Although both of these approaches seem plausible, many variations within 
these two methods (especially for the former) were tried without success to 
achieve stable results. These variations are described in Appendix B; a total 
of seven approaches were examined. Numerical instabilities due to the stiffness 
of the source-dominated k and e transport equations seemed to be the cause of 
the problem. All non-linear equations display varying degrees of stiffness 
depending on the nature of the equation, the type and strength of the non- 
linearity, and the effects of coupling with other equations in the solution 
algorithm. The k and e equations are especially stiff in the sense that they 
are coupled and at high Reynolds numbers^ they become source-dominated, i.e., 
the levels of k and e are determined to a large extent by the relative magni- 
tudes of their generation (source) and destruction (sink) rates. At high 
Reynolds numbers convective and diffusive transport are comparatively minor. 
Therefore, the treatment of these source-sink terms is of vital importance in 
achieving stable results. Our efforts so far have not produced a stable method 
for treating source-dominated equations in the MacCormack time-split finite 
difference technique of TWODLE. 

3.3.2 Scheme II 

This method adopts a new formulation where the k and e equations are solved 
only once at each time step after the velocities are advanced in time and the 


t 


Basically for all flows of aerodynamic interest. 
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source terms are calculated using these updated velocities. In this formulation 
the solution algorithm has the following form 

(1) Guess or specify initial fields of all variables. 

(2) Go through the symmetric operator sequence L (At/2) L (At) L (At/2) 

y * y 

to obtain the new (n + 1) p, U and V fields. 

(3) Solve the k and e transport equations (explicitly) using the new 
(n + 1) velocity and density, and the old (n)k, e and y^ fields to 
obtain the updated (n + 1 ) k and e values- 

(4) Calculate y. from its definition in terms of c , p, k and e. 

L ]i 

(5) Go back to step 2 until steady-state is reached. 

It should be noted that this method uses a "mixed" time level formulation 
(U, V and p are at the n + 1 level, and k, e and y^ are at the n level) and 
solves the k and e equations sequentially and explicitly. The finite difference 
form of the k and e equations and the changes introduced into TWODLE for model 
implementation are outlined in Appendix C. With this approach a stable solution 
appears to be achieved, although the computations have not yet been carried to 
steady state. 

Table 3.5 presents the U-velocity profiles as predicted by the k-e, and 

the TWODLE built-in mixing length and Baldwin-Lomax algebraic models for a 

Mach 5 10° duct compression corner flow. Results were obtained using a 21 by 

21 node mesh for a solution domain length of L = 0.1m and a width of h = 0.02m. 

A uniform inlet velocity profile and turbulence intensity were assumed. A 
* 1 * ”9 

constant At of 7.5 x 10 second was employed and 1000 time steps were taken. 

Results show that the k-e model generally predicts profiles that are in close 
agreement with the Baldwin-Lomax model. The mixing length model, on the other 
hand, produces profiles that are less steep in the vicinity of the walls than 
the other models. The k, e and y^ fields as predicted by the k-e model also 


^A constant At was used to aid the stability of the k-e model. 

j.1 

TT Computing budget considerations limited the total time steps taken to 1000. 
Thus steady-state results were not obtained. 
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TABLE 3.5. 


U-Velocity Profiles for the Mach 5 10 
Ducted Compression Corner Test Case 


x/L = 

0.20 

u/u . 

' in 



x/L = 0. 

,35 

u / u i„ 


y/h 

k-e 

ITURB = 2 

ITURB = 


y/h 

k-e 

ITUR8 = 

2 ITURB = 

0.0 

0.0 

0.0 

0.0 


0.132 

0.0 

0.0 

0.0 

0.010 

0.716 

0.699 

0.577 


0.142 

0.610 

0.606 

0.532 

0.026 

0.903 

0.901 

0.835 


0.155 

0.807 

0.806 

0.723 

0.048 

0.967 

0.967 

0.963 


0.174 

0.917 

0.916 

0.865 

0.077 

0.996 

0.997 

0.998 


0.199 

0.937 

0.937 

0.930 

0.117 

1.000 

1.000 

1.000 


0.234 

0.934 

0.961 

0.967 

0.922 

1.000 

1.000 

1.000 


0.279 

0.990 

0.991 

0.994 

0.952 

1.000 

1.000 

0.989 


0.336 

1.000 

0.999 

1.000 

0.974 

0.986 

0.985 

0.871 


0.404 

1.000 

1.000 

1.000 

0.989 

0.884 

0.874 

0.600 


0.482 

1.000 

1.000 

1.000 

1.000 

0.0 

0.0 

0.0 


0.933 

1.000 

1.000 

0.998 






0.959 

0.999 

0.999 

0.971 






0.977 

0.982 

0.981 

0.816 






0.991 

0.870 

0.866 

0.560 






1.000 

0.0 

0.0 

0.0 

Where 

h is the height of the duct, 

L is the length of the solution domain 

• afld U in 

is the maximum 

inlet 

velocity. 

ITURB is a code parameter that 

controls selection of the 

turbulence model. 

ITURB 

= 1 is a mixing length model 

for non- 

•recirculating flows 

and ITURB 

= 2 is the 

i Baldwin-Lomax 


algebraic eddy viscosity model. 


U-Velocity Profiles for the Mach 5 10° 

Ducted Compression Corner Test Case (continued) 


TABLE 3.5. 


x/L = 0 

50 

u/u. 

' in 


x/L = 

0.80 

u/u. 

' in 


y/h 

k-e 

ITURB = 2 

ITURB = 1 

y/h 

k-e 

ITURB = 2 

ITURB = 

0.176 

0.0 

0.0 

0.0 

0.176 

0.0 

0.0 

0.0 

0.185 

0.656 

0.631 

0.490 

0.185 

0.762 

0.818 

0.540 

0.198 

0.845 

0.843 

0.701 

0.198 

0.965 

0.970 

0.792 

0.216 

0.946 

0.951 

0.891 

0.216 

0.997 

0.998 

0.967 

0.240 

0.998 

0.998 

0.991 

0.240 

1.000 

1.000 

0.998 

0.272 

1.005 

1.005 

1.003 

0.272 

1.000 

1.000 

1.000 

0.315 

1.002 

1.002 

1.002 

0.315 

1.000 

1.000 

1.000 

0.369 

1.001 

1.001 

1.000 

0.369 

1.000 

1.000 

1.000 

0.434 

1.000 

1 .000 

1.000 

0.434 

1.000 

1.000 

1.000 

0.742 

1.000 

1.000 

1.000 

0.742 

1.000 

1.000 

1.000 

0.861 

1.000 

1.000 

1.000 

0.861 

1.000 

1.000 

1.000 

0.904 

1.000 

1.000 

1.000 

0.904 

1.000 

1.000 

1.000 

0.936 

1.000 

1.000 

0.998 

0.936 

1.000 

1.000 

0.998 

0.961 

0.996 

0.997 

0.967 

0.961 

0.997 

0.998 

0.967 

0.978 

0.967 

0.970 

0.793 

0.978 

0.965 

0.970 

0.792 

0.991 

0.818 

0.822 

0.541 

0.991 

0.762 

0.818 

0.540 

1.000 

0.0 

0.0 

0.0 

1.000 

0.0 

0.0 

0.0 


seem plausible. However, before any valid comparisons can be made between the 
models steady-state results should be obtained for a number of test cases 
representing both elliptic and parabolic flows. 

From these results, it appears that a stable way of incorporating the k-e 
model has been devised: no instabilities were encountered within 1000 time 

steps. However, in its present form this scheme has two drawbacks: firstly, it 

requires a small time step (approximately 2 orders of magnitude smaller than the 
time step used with the other models) and secondly, the k and e transport equa- 
tions are solved in the physical rather than in the transformed coordinates, 
which introduces inaccuracies for irregular geometries. The first drawback can 
be eliminated by devising an implicit solution algorithm (ADI, block-solver, 
etc.); ideally the k and e equations will be solved simultaneously (not sequen- 
tially) using a linearization technique such as Newton-Raphson. This will 
remove the small time step restriction and produce results that are more plausi- 
ble (since physically things do happen simultaneously). The second drawback 
can be removed by solving the k and e equations in the transformed coordinates 
of TWODLE. This is straightforward analytically, but in practice it would 
require extensive programming efforts and code changes. Both of these sugges- 
tions are important and should be considered in future code development efforts 
if this method is adopted. If, however, a new implicit version of TWODLE that 
does not use time-split finite differencing is available in the near future, the 
k and e transport equations may be successfully implemented into that solution 
algorithm since the problems with the initial version of TWODLE were largely due 
to the time-split finite differencing technique. 
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4. CONCLUSIONS AND FUTURE WORK 


The conclusions reached from the work outlined in this report can be 
summarized as follows. 

4.1 SUBSONIC PLANAR RECIRCULATING FLOWS 

(a) The relative performance of the models in subsonic plane backward- 
facing step flow computations is region-dependent, i.e., best predictions are 
not necessarily obtained by the same model in both the reverse flow and recovery 
regimes. 

(b) The "modified" ASM produces the best predictions in the reverse flow 
region but computes too slow a recovery rate in the relaxation regime where the 
other models show better agreement with data. 

(c) A modular approach that uses the "modified" ASM in the reverse flow 
region and the standard version in the recovery regime appears to be the most 
appropriate scheme in predicting subsonic plane backward-facing step flows. 

(d) The standard and multi -scale k-e models are not recommended for sub- 
sonic planar recirculating flow predictions. However, the concept of multi-time 
and -length scales is physically sound, and is also appealing in the sense that 
it enables separate modeling of each energy region in turbulent flows. 

4.2 SUBSONIC AXISYMMETRIC RECIRCULATING FLOWS 

Evaluation of the relative performance of the models in subsonic axi sym- 
metric backward-facing step flow calculations was not attempted since only one 
diameter ratio (2:1) was tested. Preliminary results indicate very little 
difference between the k-e, "modified" k-e, and ASM predictions; all show good 
agreement with the Chaturvedi (Ref. 43) data. The "modified" ASM predicts too 
long a reattachment length and consequently suffers from scaling problems. More 
diameter ratios need to be tested before the above observations can be 
general ized. 
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4.3 SUPERSONIC RECIRCULATING FLOWS 


Incorporation of the k-e model in TWODLE for supersonic recirculating flow 
predictions was hindered by numerical stability problems. Finally, a scheme was 
devised that appears to be stable but requires a very small time step due to the 
explicit nature of the scheme. Also, the fact that the k and e transport equa- 
tions are solved in the physical coordinates may introduce inaccuracies for 
irregular geometries. Two options are available: modify this scheme so that 

the k and e equations are solved implicitly and simultaneously in the trans- 
formed coordinates, or implement the k-e model into the new implicit version of 
TWODLE that does not use time-split finite differencing. 

4.4 FUTURE WORK 

Future work in the examination of turbulence models for application to 
SCRAMJET combustors is expected to concentrate in three major areas. The exami- 
nation of turbulence model performance in axisymmetric recirculating flows will 
be extended to a greater range of diameter ratios to examine whether the model 
performance observed in the work discussed in this report shows similar charac- 
teristics at other area ratios. Implementation of turbulence models in the 
supersonic-flow TWODLE code will continue, with further refinement of the k-e 
two-equation formulation and incorporation of an ASM formulation as well. 
Questions that have arisen from this work regarding the use of an explicit 
formulation for the turbulence model have ramifications not only with respect 
to computational techniques but also with respect to the architecture of the 
computer utilized, and these will be further investigated. Finally, the turbul- 
ence model assessment work will be extended to reacting flows with an emphasis 
on the development and implementation of techniques to account for turbulence- 
finite rate chemical kinetics interaction. 

All of this work has as its goal the definition of a set of turbulence 
models suitable for the different characteristic regions encountered in a scram- 
jet flowfield. Along with this definition, algorithms to allow transition from 
one turbulence model to another within the context of a given flowfield solution 
procedure are to be developed. Major strides toward the definition of suitable 
turbulence models for different flowfield regions have been made in the work 
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described in this report, and the work currently underway is expected to complete 
the turbulence model definition effort. 
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APPENDIX A 

THE k-e MODEL FORMULATION FOR THE TWODLE CODE 


The mean flow and turbulence model equations for the k-e model can be 
written as follows for two-dimensional elliptic planar flows: 

Conservation of Mass 
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where 


U = streamwise mean velocity component 
V = transverse mean velocity 
P = pressure 

k = turbulence kinetic energy 
e = turbulence kinetic energy dissipation rate 
p = density 
u = dynamic viscosity 
p = production rate of 

- 2/3k (i + 1) 

\ij = + y (total viscosity) 

M t = c y p (turbulent viscosity) 

a, and a are the Prandtl numbers for k and e, respectively, and c , c 
k e e i e 2 

and c are constants, 
u 



Following Ref. 46 these equations can be put in TWODLE form by defining the 0 , 
0, 0 and 0 vectors as 
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^ = V v 
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(%f + ($]*(%*$ -$*}-*'** 


°k ' 3/2 f 

K C k 


1 


(A10) 
(An) 
(AT 2) 


a = 


1/2 


D = 


(c - c )c 
v e 2 e ] y 


9U 3V 

-j^r + -gy compressible flows 


0 incompressible flows 

Currently recommended values for the constants are 
c^ = 0.09 

c. = 0.22 
K 1 


'e ] =1.44 


1.92 


K = 0.4187 

with these values a k and become 
°k + “ 


(A13) 


a = 1.217 

£ 


actually becomes 0.614, however a value of 1 is more common in the litera- 
ture. There are no significant differences in the predictions obtained with 
these two values. 
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APPENDIX B 


UNSUCCESSFUL APPROACHES TO THE INCORPORATION OF 
THE TWO-EQUATION k-e MODEL IN THE TWODLE CODE 

Several variations within the confines of one approach to the incorporation 
of the two-equation k-e turbulence model in the TWODLE code were examined during 
this program. Although all were ultimately unsuccessful, they are outlined in 
this Appendix so that similar problems with other stiff equation subsets may be 
avoided in numerical algorithms similar to that used in TWODLE. Testing was 
carried out for a Mach 5 10° ducted expansion-compression case using a 21 by 21 
mesh. Numerical stability problems were encountered due to the stiffness of 
the source - dominated k and e transport equations (equation (A1 ) in Appendix 
A). 

Within the TWODLE format there are at least two ways of incorporating the 

k and e equations into the solution procedure. One is to solve them sequentially 

in the l. y (At/2) L x (At) L y (At/2) symmetric operator sequence at both the interme- 
diate (n + 1) and new (n + 1) time levels with the generation and destruction 
terms calculated only during the L operator step. The second way is to solve 

A 

them only once at each time step (n to n + 1 ) after the velocities are advanced 
ih time and source terms are calculated using these updated velocities. Within 
the first mentioned scheme, seven different variations were tested to overcome 
the stiffness problem. These are discussed next. 

The seven schemes tested cover a range of source-sink treatments from fully- 

explicit to variations within a quasi-impl icit formulation. The diffusion and 

convection terms are calculated explicitly in each case. The general form of 
the source-sink terms for the k and e equations are recalled from equation (A1 ) 
as 


and 


p(p - e) for k. 


(Bl) 


pe/k(c p - c ) f° r e» 

£ I E 


(B2 ) 
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where all the terms are defined in Appendix A. Details of each treatment are 

summarized below in connection with these two expressions. In all cases a con- 

_g 

stant At of 7.5 x 10 sec was used throughout the computations. 

Case 0 

This is the fully explicit case where p, p, e, and k are all evaluated at 
the old time level . 

Calculations became unstable after 59 time steps. 

Case 1 

This is the first variation within the quasi-impl icit formulation, p is 
calculated implicitly after the new velocity field is obtained, p, c, and k are 
explicit at the old time level. 

Calculations became unstable after 59 time steps. 

Case 2 

This is the second variation within the quasi -implicit formulation, p and 

p are calculated implicitly after the new velocity field is obtained, e and k 

are explicit at the old time level. 

Calculations became unstable after 59 time steps. 

Case 3 

This is the third variation within the quasi-impl icit formulation, p and 
p are calculated implicitly after the new velocity field is obtained. The k 
that appears in equation (82) is now substituted for by the new k value, e is 
explicit at the old time level. 

Calculations became unstable after 62 time steps. 

Case 4 

This is the fourth variation within the quasi-impl icit formulation, p and 

p are calculated implicitly after the new velocity field is obtained. The term 

e/k that appears in equation (B2) is . now calculated as 

2 

c pk / y f c pk 

~ / L, - M _ M 
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where both p and k are at the new time level, e is explicit at the old time. 
Calculations became unstable after 60 time steps. 


Case 5 

This is the fifth variation within the quasi - imp! icit formulation, p and p 
are calculated implicitly after the new velocity field is obtained. The terms 
e/k and e that appear in equation (B2) are now calculated as 


e/k 


c pk /y. c pk 

U t IT 


and 


e 


c pk 2 

u 


respectively, where both p and k are at the new time level, e in equation (Bl) 
is explicit at the old time. 

Calculations become unstable after 59 time steps. 

Case 6 


This is the sixth and final variation within the quasi-impl icit formula- 
tion. p and p are calculated implicitly after the new velocity field is 
obtained. The term e in equation (Bl ) is calculated as 

e = c^pk 2 /u t 


where both p and k are at the old time level. The e/k and e terms that appear 
in equation (B2) are calculated as 


e/k 



and 


y k2 

^t 
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respectively, where both p and k are now at the new time level, 
go unstable after 27 time steps. 


Calculations 



APPENDIX C 


FINITE DIFFERENCE FORM OF THE k AND e 
TRANSPORT EQUATIONS FOR THE TWODLE CODE 

The finite difference forms given below use a "mixed" time level formulation 
where U, V, p are at the new n + 1 time level, and k, e, u t are at the old time 
level n. 



a. Temporal Term 



3t At 

b. Convective Terms (with upwind differencing) 
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if V > 0 
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if V < 0 
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Source-Sink Terms 
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d. Diffusion Terms (with central differencing) 
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where 
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c. Source-Sink Terms 
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d. 


Diffusion Terms (with central differencing) 



C.3 Wall Function Treatment 

Wall functions for the modified TWODLE code are obtained using the non- 
equilibrium wall -functions of Chieng and Launder (Ref. 26) discussed in 
Section 2.9. However, since TWODLE uses nodal values rather than control 
volume averages 1 - , some changes were made in implementing these wall functions. 

A typical near-wall region is shown in Figure C.l. Here it is assumed that node 
w is at the wall, and w + 1 is in the fully turbulent region 


w + 


l k v 


1/2 


20 , 


where y, , , is the distance from node w + 1 to the wall, and k is the turbu- 
w * i v 

lent kinetic energy at the edge of the viscous-sublayer, y . Following 

Ref. 26 the wall shear stress T can be expressed as 

w 


T w = 


K*pU 


w + 




InE* 


r V4 + 


l k v 


1/2 


Chieng and Launder wal 1 -functions use near-wall cell integration to calculate 
mean production and dissipation rates for those cells. 
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1/2 K*Re , 

where k* = 0.4187c^ and E* = e /Re y . Also the shear stress at node 
w + 2 by definition is 


T. 


. /9U , 9V \ 

w + 2-May 9x/ w + 2 


If it can be conjectured that the shear stress varies linearly between the wall 
and node w + 2 (a plausible assumption for most flows with a reasonably fine 
near-wall grid), the shear stress at node w + 1 can be obtained by interpolating 
between nodes w and w + 2 


T w + 1 T w 


+ ( T w + 2 - T w ) 


\ + 1 - yj 


(y« + 2 - *«) 

Then the turbulent viscosity at w + 1 can be calculated from its definition 
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The production rate of turbulent kinetic energy at w + 1 now becomes 
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and the near-wall dissipation rates are expressed as 
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/ 3U 9V_\ 
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and 




following Spalding (Ref. 29), and Pope and Whitelaw (Ref. 28), respectively. 


'The universal viscous-sublayer thickness constant Re y is assumed to be 20. 
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